Method for calculating the visible window of an optical satellite for ground targets

Abstract

This invention relates to a method for calculating the visibility window of an optical satellite over ground targets, comprising the following steps: considering the orbital attitude motion of the optical satellite and the installation and field-of-view parameters of the onboard optical camera, establishing a visibility determination function for ground point targets by the optical satellite; fitting the determination function with a piecewise cubic polynomial function to calculate the time window for visibility of ground point targets by the optical satellite; using a discretization method to calculate the time window for visibility of ground line targets by the optical satellite; considering the motion law of the optical satellite, calculating the time window for visibility of ground triangular region targets by the optical satellite; and using a polygon triangulation method to calculate the time window for visibility of ground targets in any polygonal region.

Description

Technical Field

[0001] This invention relates to the field of optical satellite technology, and in particular to a method for calculating the visible time window of an optical satellite for ground targets. Background Technology

[0002] With the development of aerospace technology, optical satellite-based Earth remote sensing observation technology has been widely applied in many fields such as meteorological observation, environmental monitoring, disaster early warning, and urban planning. Rapidly and accurately calculating the visible time window of optical satellites over ground targets is fundamental to scheduling optical remote sensing satellite observation missions. Ground targets can be point targets such as airports and ports, line targets such as rivers and borders, or arbitrary polygonal areas. Furthermore, during the operation of the optical camera, the optical satellite may also perform attitude maneuvers to achieve complex modes of Earth remote sensing observation. Meeting these practical application requirements and rapidly and accurately calculating the visible time window of optical satellites over ground targets is a relatively complex problem.

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

[0004] To address the technical problems existing in the prior art, the present invention aims to provide a method for calculating the visible window of ground targets by an optical satellite, which can realize the rapid and accurate calculation of ground point, line, and area targets by an optical satellite.

[0005] To achieve the above-mentioned objective, this invention provides a method for calculating the visibility window of an optical satellite over 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, establish a function for determining the visibility of the optical satellite to the ground point target. Fit the function with a piecewise cubic polynomial function and 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 during which the optical satellite is visible to the ground line target, and considering the motion law of the optical satellite, calculate the time window during which the optical satellite is visible to the ground triangular region target;

[0009] Step S4: Based on the visible time window of the optical satellite for the ground triangular region target, the visible time window of the optical satellite for any polygonal region target on the ground is calculated using the polygonal triangulation method.

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

[0011] Step S11: Based on the optical satellite's Earth-fixed system ephemeris, calculate the Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite at any time within the time interval of the optical satellite's Earth-fixed system ephemeris.

[0012] Step S12: Based on the optical satellite attitude sequence, calculate the attitude data of the optical satellite at any time within the time interval of the optical satellite's Earth-Fixed System ephemeris.

[0013] Step S13: Based on the optical satellite attitude data and the optical satellite optical camera installation parameters, calculate the transformation matrix from the Earth-fixed system to the camera installation coordinate system of the optical satellite at any given time;

[0014] Step S14: Considering the attitude motion of the optical satellite and the installation and field of view parameters of the optical camera of the optical satellite, establish a determination function for the visibility of the ground point target by the optical satellite;

[0015] Step S15: Establish the first and second derivatives of the decision function based on the central difference quotient;

[0016] Step S16: Using an adaptive step size strategy, construct a piecewise cubic Hermitian interpolation polynomial function to fit the decision function;

[0017] Step S17: Find the roots of the piecewise cubic Hermitian interpolation polynomial function and calculate the time window set that satisfies the decision function;

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

[0019] According to one technical solution of the present invention, in step S11, the optical satellite Earth-Fixed System ephemeris is represented as follows: , This indicates the ephemeris time point in the Earth-Fixed System ephemeris of the optical satellite. and This indicates the ephemeris time node in the Earth-Fixed System ephemeris of the optical satellite. The corresponding Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite;

[0020] For any time The Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite are expressed as follows: and The value is obtained by segmenting five times using Lagrange interpolation; The corresponding ephemeris time node in the optical satellite Earth-Fixed System is: Let k represent the degree of Lagrange interpolation, then we have

[0021]

[0022] in, Let be the basis function of the k-th Lagrange interpolation, expressed as:

[0023] ;

[0024] In step S12, the optical satellite attitude sequence is represented as follows: ,in for The time-based optical satellite body coordinate system is a unit quaternion relative to the VVLH coordinate system. The first three are the quaternion vector parts, and the fourth is the quaternion scalar part.

[0025] For any time The attitude data of the optical satellite is represented as follows: The attitude data of the optical satellite is obtained through piecewise spherical linear interpolation based on the optical satellite attitude sequence. The time node in the optical satellite attitude sequence corresponding to the time is Then there is

[0026]

[0027] Where k0, k1, and θ are the initial weight, the final weight, and the angle between the two quaternions, respectively, denoted as:

[0028] ;

[0029] In step S13, the transformation matrix from the Earth-fixed coordinate system to the camera mounting coordinate system of the optical satellite is expressed as follows:

[0030]

[0031] in, It is the camera mounting matrix, which is a constant matrix; This is the transformation matrix from the VVLH coordinate system to the optical satellite body coordinate system. This is the transformation matrix from the Earth-fixed system to the VVLH coordinate system of the optical satellite. Depend on Attitude data of the optical satellite at the specified time The calculation is obtained and expressed as

[0032]

[0033]

[0034] The Earth-fixed position vector of the optical satellite at time t and the velocity vector of the ground-solid system The calculation is obtained and expressed as

[0035]

[0036]

[0037]

[0038]

[0039] in, , , and These represent the coordinates of the unit vectors of the VVLH coordinate axes in the Earth-fixed system. , This is the Earth's rotational angular velocity vector.

[0040] According to one technical solution of the present invention, step S14 specifically includes:

[0041] Step S141: Obtain the geodetic coordinates of the ground point target. The ground point target is then converted into rectangular coordinates to obtain the ground-fixed system vector of the ground point target. , represented as

[0042]

[0043] in, Longitude Latitude of the earth denoted as the ground height; x, y, and z represent the rectangular coordinates of the ground-fixed vector of the target point. The radius of curvature of the circle is denoted as . , These are the semi-major axis and eccentricity of the Earth's ellipsoid, respectively.

[0044] Step S142: Transform the ground-fixed coordinate system vector of the ground point target to the camera mounting coordinate system, and construct the determination function of the ground point target within the field of view of the optical satellite, expressed as:

[0045]

[0046]

[0047] in, and These 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) represent the coordinate components of the vector from the optical satellite to the ground point target in the camera mounting coordinate system; and The function for determining whether a ground point target is within the field of view of the optical satellite;

[0048] Step S143: Construct the minimum elevation angle determination function for the optical satellite relative to a ground point target, expressed as:

[0049]

[0050] Where, n p This represents the geodetic horizontal plane normal vector of the ground point target. ; The minimum elevation angle of the optical satellite to the ground point target;

[0051] Step S144: Construct the ground solar altitude angle determination function, expressed as:

[0052]

[0053] in, This is the minimum solar altitude angle at ground level; Let be the Earth-fixed system solar position vector at time t;

[0054] Step S145: Simultaneously satisfy the judgment function At that time, the optical satellite becomes visible to the ground point target.

[0055] According to a technical solution of the present invention, in step S15, the first and second derivatives of the function for determining the visibility of ground point targets by an optical satellite based on the central difference quotient are expressed as follows:

[0056]

[0057] in, The step size for the difference quotient is fixed at 0.1 seconds.

[0058] According to one technical solution of the present invention, in step S16, constructing the piecewise cubic Hermitian interpolation polynomial function specifically includes:

[0059] Step S1601: Set the maximum interpolation step size It is one-quarter of the orbital period;

[0060] Step S1602: Set the current interpolation node Calculate the decision function in Function value and derivative value at point , , Save the current interpolation node data , , , ;

[0061] Step S1603: Calculate the interpolation step size ;

[0062] Step S1604: Set the successor interpolation node Calculate the decision function in Function value and derivative value at point , , ;

[0063] Step S1605: In the interpolation interval Construct a two-point cubic Hermitian interpolation polynomial. ,calculate The second derivative value at the endpoints of the interpolation interval , ;

[0064] Step S1606: Calculate the fitting error ,like Then shorten the interpolation step size Proceed to step S1604;

[0065] Step S1607, Calculation In the open interval If the equation has roots, then the equation has roots. and Then adjust the interpolation step size. Proceed to step S1604;

[0066] Step S1608: Save the subsequent interpolation node data , , , ;

[0067] Step S1609: Update the current interpolation node data. , , , ;

[0068] Step S1610, if If the result is positive, proceed to step S1603; otherwise, end the calculation process and obtain the difference node set. , where n represents the number of interpolation nodes.

[0069] According to one technical solution of the present invention, step S17 specifically includes:

[0070] Step S171: Set the time window flag isInWin=false;

[0071] Step S172: Set the difference node set The first interpolation node data in the array is used as the current interpolation node data. , , Set window start time Set isInWin=true;

[0072] Step S173: Obtain subsequent interpolation node data , ;

[0073] Step S174, if ,calculate exist inner root Otherwise, proceed to Step 5; if isInWin is true, set the window end time. It also saves the window data; if isInWin is false, it sets the window start time. Modify the time window flag: isInWin = !isInWin;

[0074] Step S175: Update the current interpolation node data. , ;

[0075] Step S176, if If yes, proceed to step S173; otherwise, proceed to step S177.

[0076] Step S177: If isInWin is true, then set the window end time. And save the window data.

[0077] According to a technical solution of the present invention, in step S2, the visible window algorithm for ground line targets by the optical satellite includes:

[0078] Let the two endpoints of the ground line target be respectively , The ground line targets are discretized into a set of ground point targets by using an equidistant partitioning method along the geodesic line. The algorithm for the visibility window of ground point targets by optical satellites is applied to obtain the visibility time window of each ground point target in the set of ground point targets by the optical satellites. Then, a union calculation is performed to obtain the visibility time window of the ground line targets by the optical satellites, expressed as:

[0079]

[0080] .

[0081] According to a technical solution of the present invention, in step S3, the optical satellite's visible window algorithm for the ground triangular region specifically includes:

[0082] Applying the optical satellite's visibility window algorithm for ground line targets to the three sides of the aforementioned ground triangular region target, three time window sets can be obtained, represented as follows:

[0083]

[0084] First, perform a union operation on these three time window sets, which is represented as:

[0085]

[0086] Then iterate through For each time window in the data, for two adjacent time windows... , ,like If the time window is less than half an orbital period, then two adjacent time windows will be merged into one time window. ;

[0087] Repeat the above process until... Within two adjacent time windows, there is no longer any In cases where the time window is less than half an orbital period, the final time window is obtained. This refers to the time window during which the optical satellite can see the target in the triangular area on the ground.

[0088] According to one technical solution of the present invention, step S4 specifically includes:

[0089] Let G be an arbitrary polygonal region on the ground. ,in The coordinates of the polygon vertices are the geodetic coordinates.

[0090] Using the latitude and longitude of the vertices of the polygonal region as the horizontal and vertical coordinates respectively, the polygonal region is divided into sections using a planar polygon triangulation method. A triangular region, represented as

[0091]

[0092] in , , The coordinates of the vertices of the triangular region are given.

[0093] By applying the aforementioned optical satellite visibility window algorithm for triangular regions on the ground, a set of time windows for each triangular region to be visible to the optical satellite can be obtained. This set is then combined to obtain the time window for the optical satellite to be visible to targets in any polygonal region on the ground, denoted as...

[0094]

[0095] .

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

[0097] The present invention provides a method for calculating the visibility window of an optical satellite over ground targets. First, considering the orbital attitude motion of the optical satellite and the installation and field-of-view parameters of its optical camera, a visibility determination function for ground point targets is established. This function is then fitted using a piecewise cubic polynomial function to calculate the visibility time window for ground point targets. Second, based on the visibility window algorithm for ground point targets, a discretization method is used to calculate the visibility time window for ground line targets. Third, based on the visibility window algorithm for ground line targets, considering the motion characteristics of the optical satellite, the visibility time window for ground triangular region targets is calculated. Finally, based on the visibility window algorithm for ground triangular regions, a polygonal triangulation method is used to calculate the visibility time window for ground targets in any polygonal region. This invention considers the attitude maneuvering of the optical satellite and the pointing constraints of the payload installation, employing a unified technical framework to achieve rapid and accurate calculation of ground point, line, and region targets by the optical satellite. This solves the problem that the calculation accuracy and efficiency of traditional tracking and propagation methods heavily rely on the selected step size and have low computational efficiency. Attached Figure Description

[0098] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly described below. Obviously, the drawings described below are merely some embodiments of the present invention, and those skilled in the art can obtain other drawings based on these drawings without creative effort.

[0099] Figure 1 This schematic diagram illustrates a flowchart of a method for calculating the visibility window of an optical satellite for ground targets according to an embodiment of the present invention.

[0100] Figure 2 This schematic diagram illustrates the specific flowchart for calculating the visible time window of an optical satellite for a ground point target according to an embodiment of the present invention.

[0101] Figure 3 This is a schematic diagram illustrating the principle of optical satellite visibility determination of ground point targets according to an embodiment of the present invention.

[0102] Figure 4 This schematic diagram illustrates the construction flowchart of a piecewise cubic interpolation polynomial provided in an embodiment of the present invention.

[0103] Figure 5 This schematic diagram illustrates the flowchart for generating the time window for an optical satellite to be visible to a ground point target, provided in an embodiment of the present invention.

[0104] Figure 6 This schematic diagram illustrates the discretization of ground line targets into a set of point targets according to an embodiment of the present invention.

[0105] Figure 7 This is a schematic diagram illustrating the division of an arbitrary polygon on the ground into triangles according to an embodiment of the present invention. Detailed Implementation

[0106] The description of the embodiments in this specification should be taken in conjunction with the accompanying drawings, which should form part of the complete specification. In the drawings, the shape or thickness of the embodiments may be exaggerated and may be indicated in a simplified or convenient manner. Furthermore, parts of the various structures in the drawings will be described separately; it is worth noting that elements not shown in the figures or not described in words are in a form known to those skilled in the art.

[0107] The descriptions of the embodiments herein, including any references to directions and orientations, are for ease of description only and should not be construed as limiting the scope of the invention. The following description of preferred embodiments involves combinations of features, which may exist independently or in combination; the invention is not particularly limited to the preferred embodiments. The scope of the invention is defined by the claims.

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

[0109] 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, establish the visibility determination function of the optical satellite to ground point targets, use a piecewise cubic polynomial function to fit the determination function, and calculate the time window when the optical satellite is visible to ground point targets.

[0110] Step S1 specifically includes:

[0111] Step S11: Based on the optical satellite's Earth-fixed system ephemeris, calculate the optical satellite's Earth-fixed system position vector and Earth-fixed system velocity vector at any time within the time interval of the optical satellite's Earth-fixed system ephemeris.

[0112] This embodiment does not impose any constraints on the satellite orbit type, nor does it involve the calculation of the satellite orbit. Instead, it directly uses the satellite's Earth-Fixed System ephemeris to describe the satellite's orbital motion state.

[0113] In step S11, the satellite ephemeris is within a time interval. The set of discrete satellite position and velocity vectors within the Earth-Fixed System is represented by the optical satellite's Earth-Fixed System ephemeris as follows: , This indicates the ephemeris time point in the Earth-Fixed System ephemeris of an optical satellite. and Indicates the ephemeris time nodes in the ephemeris of the Earth-Fixed System of optical satellites. The corresponding Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite;

[0114] For any time The Earth-fixed system position vector and Earth-fixed system velocity vector of an optical satellite are expressed as follows: and The value is obtained by segmenting five times using Lagrange interpolation; The corresponding ephemeris time node in the optical satellite Earth-Fixed System is Let k represent the degree of Lagrange interpolation, then we have

[0115]

[0116] in, Let be the basis function of the k-th Lagrange interpolation, expressed as:

[0117] ;

[0118] Step S12: Based on the optical satellite attitude sequence, calculate the attitude data of the optical satellite at any time within the time interval of the Earth-Fixed System ephemeris of the optical satellite.

[0119] This embodiment does not impose any constraints on the satellite attitude mode, nor does it involve the calculation of the satellite attitude. Instead, it directly uses the satellite's own system attitude sequence relative to the VVLH coordinate system to describe the satellite's attitude motion state.

[0120] In step S12, the satellite attitude sequence is within a time interval An internal set of discrete satellite attitude quaternions, the optical satellite attitude sequence is represented as: ,in for The time-based optical satellite body coordinate system is a unit quaternion relative to the VVLH coordinate system. The first three are the quaternion vector parts, and the fourth is the quaternion scalar part.

[0121] For any time The attitude data of the optical satellite is represented as The attitude data of the optical satellite is obtained through piecewise spherical linear interpolation based on the optical satellite attitude sequence. The time node in the optical satellite attitude sequence corresponding to the time is Then there is

[0122]

[0123] Where k0, k1, and θ are the initial weight, the final weight, and the angle between the two quaternions, respectively, denoted as:

[0124] ;

[0125] Step S13: Based on the optical satellite attitude data and the optical camera installation parameters, calculate the transformation matrix from the Earth-fixed coordinate system to the camera installation coordinate system of the optical satellite at any given time.

[0126] In step S13, the transformation matrix from the Earth-fixed coordinate system to the camera mounting coordinate system of the optical satellite is expressed as follows:

[0127]

[0128] in, It is the camera mounting matrix, which is a constant matrix; This is the transformation matrix from the VVLH coordinate system to the optical satellite body coordinate system. This is the transformation matrix from the Earth-fixed coordinate system to the VVLH coordinate system for optical satellites. Depend on Attitude data of the time-of-flight optical satellite The calculation is obtained and expressed as

[0129]

[0130]

[0131] The Earth-fixed system position vector of the optical satellite at time t and the velocity vector of the ground-solid system The calculation is obtained and expressed as

[0132]

[0133]

[0134]

[0135]

[0136] in, , and These represent the coordinates of the unit vectors of the VVLH coordinate axes in the Earth-fixed system. This represents the sum of the satellite's velocity vector in the Earth-solid system and the tangential velocity vector of the Earth's rotation. , This is the Earth's rotational angular velocity vector.

[0137] Step S14: Considering the attitude motion of the optical satellite and the installation and field of view parameters of the optical camera, establish a function to determine whether the optical satellite is visible to ground point targets;

[0138] Step S14 specifically includes:

[0139] Step S141: Obtain the geodetic coordinates of the ground point target. The ground point target's geodetic coordinates are converted to rectangular coordinates to obtain the ground-fixed system vector of the ground point target. , represented as

[0140]

[0141] in, Longitude Latitude of the earth y represents the ground height; x, y, and z represent the rectangular coordinates of the ground point target's ground-fixed vector. The radius of curvature of the circle is denoted as . , These are the semi-major axis and eccentricity of the Earth's ellipsoid, respectively.

[0142] Step S142: Transform the ground-fixed coordinate system vector of the ground point target to the camera mounting coordinate system, and construct the determination function of the ground point target within the optical satellite's field of view, expressed as:

[0143]

[0144]

[0145] in, and These represent the vertical half-angle and horizontal half-angle of the field of view of an optical satellite, respectively; r y (t), r x (t), r z (t) represent the coordinate components of the vector from the optical satellite to the ground point target in the camera mounting coordinate system; and This is the function for determining whether a ground point target is within the field of view of an optical satellite.

[0146] Step S143: Construct the minimum elevation angle determination function for the optical satellite relative to a ground point target, expressed as:

[0147]

[0148] Where, n p The normal vector of the ground point target to the horizontal plane. ; This is the minimum elevation angle of the optical satellite to a ground point target;

[0149] Step S145: Construct the ground solar altitude angle determination function, expressed as:

[0150]

[0151] in, This is the minimum solar altitude angle at ground level; Let be the Earth-fixed system solar position vector at time t;

[0152] Step S145: Simultaneously satisfy the judgment function At that time, optical satellites can see ground point targets.

[0153] Step S15: Establish the first and second derivatives of the decision function based on the central difference quotient;

[0154] In step S15, the first and second derivatives of the visibility determination function for ground point targets by the optical satellite are established based on the central difference quotient, expressed as follows:

[0155]

[0156] in, The step size for the difference quotient is fixed at 0.1 seconds.

[0157] Step S16: Using an adaptive step size strategy, construct a piecewise cubic Hermitian interpolation polynomial function fitting decision function;

[0158] Decision function Its form is complex, and its complete analytical expression is difficult to determine within a time interval. For the decision equation Finding roots can only be done using mechanical search methods, which involve a large amount of computation and make it difficult to guarantee that no roots are missed.

[0159] This embodiment uses two-point cubic Hermitian interpolation to construct a piecewise cubic polynomial. Function fitting decision function. Construct interpolation node set. , making Conditions met:

[0160] a.

[0161] b.

[0162] c.

[0163] Therefore, the equation is determined. The problem of finding the roots can be transformed into the problem of finding the roots of interpolation polynomials, which can effectively reduce the difficulty of solving the problem and improve the efficiency of solving the problem.

[0164] In the time interval Select a reasonable number of interpolation points to make the interpolation function It can accurately approximate the decision function To maximize computational efficiency, this embodiment constructs interpolation nodes using the following principles:

[0165] 1) Use the difference of the second derivatives at the endpoints of the interpolation interval to control the interpolation step size.

[0166] In the interval The fitting error is defined above and expressed as:

[0167]

[0168]

[0169] Given error limit ,like This shortens the interpolation interval. .

[0170] 2) Ensure that within a single interpolation interval It is a monotonic function

[0171] Given interpolation interval Construct the interpolation polynomial, expressed as

[0172]

[0173] right Differentiate, and we get

[0174]

[0175] Find the equation of the quadratic polynomial. In the open interval The roots within the interval. If the equation has no roots, then the interpolating polynomial lies within the interval. If the equation is monotonic, it accepts the interpolation interval; if the equation has roots... Then As the interpolation interval.

[0176] Following the aforementioned principles, such as Figure 4 As shown, in step S16, constructing the piecewise cubic Hermitian interpolation polynomial function specifically includes:

[0177] Step S1601: Set the maximum interpolation step size It is one-quarter of the orbital period;

[0178] Step S1602: Set the current interpolation node Calculate the decision function in Function value and derivative value at point , , Save the current interpolation node data , , , ;

[0179] Step S1603: Calculate the interpolation step size ;

[0180] Step S1604: Set the successor interpolation node Calculate the decision function in Function value and derivative value at point , , ;

[0181] Step S1605: In the interpolation interval Construct a two-point cubic Hermitian interpolation polynomial. ,calculate The second derivative value at the endpoints of the interpolation interval , ;

[0182] Step S1606: Calculate the fitting error ,like Then shorten the interpolation step size Proceed to step S1604;

[0183] Step S1607, Calculation In the open interval If the equation has roots, then the equation has roots. and Then adjust the interpolation step size. Proceed to step S1604;

[0184] Step S1608: Save the subsequent interpolation node data , , , ;

[0185] Step S1609: Update the current interpolation node data. , , , ;

[0186] Step S1610, if If the result is positive, proceed to step S1603; otherwise, end the calculation process and obtain the difference node set. , where n represents the number of interpolation nodes.

[0187] Step S17: Find the roots of the piecewise cubic Hermitian interpolation polynomial function and calculate the time window set that satisfies the criterion function;

[0188] In each interpolation interval inside, if If there is a root in this interval, then there is exactly one root, which can be found using the cubic polynomial root-finding formula. By traversing all interpolation nodes, the following can be calculated: exist Generate a set of time windows that satisfy the decision function from all roots within the set.

[0189] Step S17 specifically includes:

[0190] Step S171: Set the time window flag isInWin=false;

[0191] Step S172: Set the difference node set The first interpolation node data in the array is used as the current interpolation node data. , , Set window start time Set isInWin=true;

[0192] Step S173: Obtain subsequent interpolation node data , ;

[0193] Step S174, if ,calculate exist inner root Otherwise, proceed to Step 5; if isInWin is true, set the window end time. It also saves the window data; if isInWin is false, it sets the window start time. Modify the time window flag: isInWin = !isInWin;

[0194] Step S175: Update the current interpolation node data. , ;

[0195] Step S176, if If yes, proceed to step S173; otherwise, proceed to step S177.

[0196] Step S177: If isInWin is true, then set the window end time. And save the window data.

[0197] Step S18: Perform set operations on the set that satisfies the decision function to obtain the time window in which the optical satellite is visible to ground point targets.

[0198] By fitting four piecewise cubic polynomial functions to determine the visibility of ground point targets by optical satellites, four time window sets can be obtained.

[0199]

[0200] By performing an intersection operation on these four time window sets, the visible time window for ground point targets by optical satellites can be obtained.

[0201]

[0202] Step S2: Based on the time window in which optical satellites are visible to ground point targets, the discretization method is used to calculate the time window in which optical satellites are visible to ground line targets.

[0203] In step S2, the optical satellite's visible window algorithm for ground line targets includes:

[0204] Let the two endpoints of the ground line target be respectively , The ground line targets are discretized into a set of ground point targets by using an equidistant partitioning method along the geodesic line. By applying the optical satellite visibility window algorithm for ground point targets, the visibility time window for each ground point target in the set of ground point targets is obtained, and the union of these windows is calculated to obtain the visibility time window for ground line targets by the optical satellite, which is expressed as:

[0205]

[0206] .

[0207] Step S3: Based on the time window when the optical satellite is visible to ground targets, and considering the motion law of the optical satellite, calculate the time window when the optical satellite is visible to targets in the triangular area of ​​the ground.

[0208] In step S3, the optical satellite's visible window algorithm for the ground triangular region specifically includes:

[0209] Applying the optical satellite visibility window algorithm for ground line targets to the three sides of the ground triangular region target, three time window sets can be obtained, represented as follows:

[0210]

[0211] First, perform a union operation on these three time window sets, which is represented as:

[0212]

[0213] Then iterate through For each time window in the data, for two adjacent time windows... , ,like If the time window is less than half an orbital period, then two adjacent time windows will be merged into one time window. ;

[0214] Repeat the above process until... Within two adjacent time windows, there is no longer any In cases where the time window is less than half an orbital period, the final time window is obtained. It refers to the time window during which an optical satellite can see a target in a triangular area on the ground.

[0215] Step S4: Based on the time window in which the optical satellite can see targets in a triangular region on the ground, the polygon triangulation method is used to calculate the time window in which the optical satellite can see targets in any polygonal region on the ground.

[0216] Step S4 specifically includes:

[0217] Let G be an arbitrary polygonal region on the ground. ,in The coordinates of the polygon vertices are the geodetic coordinates.

[0218] Using the latitude and longitude of the vertices of the polygonal region as the horizontal and vertical coordinates respectively, the polygonal region is divided into sections using the planar polygon triangulation method. A triangular region, represented as

[0219]

[0220] in , , The coordinates of the vertices of the triangular region are given.

[0221] By applying the optical satellite visibility window algorithm for triangular regions on the ground, a set of time windows for each triangular region can be obtained. These time windows are then combined to obtain the time window for the optical satellite to be visible to targets in any polygonal region on the ground, denoted as:

[0222] .

[0223] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or terminal device that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or terminal device. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or terminal device that includes said element.

[0224] Finally, it should be noted that the above description represents a preferred embodiment of the present invention. It should be pointed out that although preferred embodiments have been described, 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 described herein. These improvements and modifications should also be considered within the scope of protection of the present invention. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the embodiments of the present invention.

Claims

1. A method for calculating optical satellite-to-ground target visibility windows, characterized in that, Includes the following steps: 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, establish a function for determining the visibility of the optical satellite to the ground point target. Fit the function with a piecewise cubic polynomial function and 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 during which the optical satellite is visible to the ground line target, and considering the motion law of the optical satellite, calculate the time window during which the optical satellite is visible to the ground triangular region target; Step S4: Based on the visible time window of the optical satellite for the ground triangular region target, the visible time window of the optical satellite for any polygonal region target on the ground is calculated using the polygonal triangulation method; Step S1 specifically includes: Step S11: Based on the optical satellite's Earth-fixed system ephemeris, calculate the Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite at any time within the time interval of the optical satellite's Earth-fixed system ephemeris. Step S12: Based on the optical satellite attitude sequence, calculate the attitude data of the optical satellite at any time within the time interval of the optical satellite's Earth-Fixed System ephemeris. Step S13: Based on the optical satellite attitude data and the optical satellite optical camera installation parameters, calculate the transformation matrix from the Earth-fixed system to the camera installation coordinate system of the optical satellite at any given time; Step S14: Considering the attitude motion of the optical satellite and the installation and field of view parameters of the optical camera of the optical satellite, establish a determination function for the visibility of the ground point target by the optical satellite; Step S15: Establish the first and second derivatives of the decision function based on the central difference quotient; Step S16: Using an adaptive step size strategy, construct a piecewise cubic Hermitian interpolation polynomial function to fit the decision function; Step S17: Find the roots of the piecewise cubic Hermitian interpolation polynomial function and calculate the time window set that satisfies the decision function; Step S18: Perform set operation on the set that satisfies the determination function to obtain the time window in which the optical satellite is visible to the ground point target; Step S14 specifically includes: Step S141, obtaining the geodetic coordinates of the ground point target and converting the ground point target into a rectangular coordinate to obtain the geodetic vector of the ground point target , expressed as in, Longitude Latitude of the earth denoted as the ground height; x, y, and z represent the rectangular coordinates of the ground-fixed vector of the target point. The radius of curvature of the circle is denoted as . , These are the semi-major axis and eccentricity of the Earth's ellipsoid, respectively. Step S142: Transform the ground-fixed coordinate system vector of the ground point target to the camera mounting coordinate system, and construct the determination function of the ground point target within the field of view of the optical satellite, expressed as: Where r(t) represents the determination function for ground point targets within the field of view of the optical satellite. This represents the transformation matrix from the Earth-fixed coordinate system to the camera mounting coordinate system of the optical satellite. This represents the Earth-fixed position vector of the optical satellite. and These 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) represent the coordinate components of the vector from the optical satellite to the ground point target in the camera mounting coordinate system; and The function for determining whether a ground point target is within the field of view of the optical satellite; Step S143: Construct the minimum elevation angle determination function for the optical satellite relative to a ground point target, expressed as: Where, n p This represents the geodetic horizontal plane normal vector of the ground point target. ; The minimum elevation angle of the optical satellite to the ground point target; Step S144: Construct the ground solar altitude angle determination function, expressed as: in, This is the minimum solar altitude angle at ground level; Let be the Earth-fixed system solar position vector at time t; Step S145: Satisfy the judgment function At that time, the optical satellite is visible to the ground point target. The step size of the difference quotient.

2. The method for calculating the visible window of an optical satellite for ground targets according to claim 1, characterized in that, In step S11, the optical satellite Earth-Fixed System ephemeris is represented as follows: , This indicates the ephemeris time point in the Earth-Fixed System ephemeris of the optical satellite. and This indicates the ephemeris time node in the Earth-Fixed System ephemeris of the optical satellite. The corresponding Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite; For any time The Earth-fixed system position vector and Earth-fixed system velocity vector of the optical satellite are expressed as follows: and The value is obtained by segmenting five times using Lagrange interpolation; The corresponding ephemeris time node in the optical satellite Earth-Fixed System is: Let k represent the degree of Lagrange interpolation, then we have in, Let be the basis function of the k-th Lagrange interpolation, expressed as: ； In step S12, the optical satellite attitude sequence is represented as follows: ,in for The time-based optical satellite body coordinate system is a unit quaternion relative to the VVLH coordinate system. The first three are the quaternion vector parts, and the fourth is the quaternion scalar part. For any time The attitude data of the optical satellite is represented as follows: The attitude data of the optical satellite is obtained through piecewise spherical linear interpolation based on the optical satellite attitude sequence. The time node in the optical satellite attitude sequence corresponding to the time is Then there is Where k0, k1, and θ are the initial weight, the final weight, and the angle between the two quaternions, respectively, denoted as: ； In step S13, the transformation matrix from the Earth-fixed coordinate system to the camera mounting coordinate system of the optical satellite is expressed as follows: in, It is the camera mounting matrix, which is a constant matrix; This is the transformation matrix from the VVLH coordinate system to the optical satellite body coordinate system. This is the transformation matrix from the Earth-fixed system to the VVLH coordinate system of the optical satellite. Depend on Attitude data of the optical satellite at the specified time The calculation yields the following result, which is expressed as: The Earth-fixed position vector of the optical satellite at time t and the velocity vector of the ground-solid system The calculation yields the following result, which is expressed as: in, , , These represent the coordinates of the unit vectors of the VVLH coordinate axes in the Earth-fixed system. This represents the sum of the satellite's velocity vector in the Earth-solid system and the tangential velocity vector of the Earth's rotation. , This is the Earth's rotational angular velocity vector.

3. The method for calculating the visible window of an optical satellite for ground targets according to claim 2, characterized in that, In step S15, the first and second derivatives of the decision function are established based on the central difference quotient, expressed as follows: in, The step size for the difference quotient is fixed at 0.1 seconds.

4. The method for calculating the visible window of an optical satellite for ground targets according to claim 3, characterized in that, In step S16, constructing the piecewise cubic Hermitian interpolation polynomial function specifically includes: Step S1601: Set the maximum interpolation step size It is one-quarter of the orbital period; Step S1602: Set the current interpolation node Calculate the decision function in Function value and derivative value at point , , Save the current interpolation node data , , , ; Step S1603: Calculate the interpolation step size ; Step S1604: Set the successor interpolation node Calculate the decision function in Function value and derivative value at point , , ; Step S1605: In the interpolation interval Construct a two-point cubic Hermitian interpolation polynomial. ,calculate The second derivative value at the endpoints of the interpolation interval , ; Step S1606: Calculate the fitting error ,like Then shorten the interpolation step size Proceed to step S1604; Indicates the error limit; Step S1607, Calculation In the open interval If the equation has roots, then the equation has roots. and Then adjust the interpolation step size. Proceed to step S1604; Step S1608: Save the subsequent interpolation node data , , , ; Step S1609: Update the current interpolation node data. , , , ; Step S1610, if If the result is positive, proceed to step S1603; otherwise, end the calculation process and obtain the difference node set. , where n represents the number of interpolation nodes.

5. The method for calculating the visible window of an optical satellite for ground targets according to claim 4, characterized in that, Step S17 specifically includes: Step S171: Set the time window flag isInWin=false; Step S172: Set the difference node set The first interpolation node data in the array is used as the current interpolation node data. , , Set window start time Set isInWin=true; Step S173: Obtain subsequent interpolation node data , ; Step S174, if ,calculate exist inner root Otherwise, proceed to Step 5; if isInWin is true, set the window end time. It also saves the window data; if isInWin is false, it sets the window start time. Modify the time window flag: isInWin = !isInWin; Step S175: Update the current interpolation node data. , ; Step S176, if If yes, proceed to step S173; otherwise, proceed to step S177. Step S177: If isInWin is true, then set the window end time. And save the window data.

6. The method for calculating the visible window of an optical satellite for ground targets according to claim 1, characterized in that, In step S2, calculating the time window in which the optical satellite is visible to ground line targets includes: Let the two endpoints of the ground line target be respectively , The ground line targets are discretized into a set of ground point targets by using an equidistant partitioning method along the geodesic line. The algorithm for the visibility window of ground point targets by optical satellites is applied to obtain the visibility time window of each ground point target in the set of ground point targets by the optical satellites. Then, a union calculation is performed to obtain the visibility time window of the ground line targets by the optical satellites, expressed as: 。 7. The method for calculating the visible window of an optical satellite for ground targets according to claim 6, characterized in that, In step S3, calculating the time window during which the optical satellite is visible to targets in the triangular region on the ground specifically includes: Applying the optical satellite's visibility window algorithm for ground line targets to the three sides of the aforementioned ground triangular region target, three time window sets can be obtained, represented as follows: First, perform a union operation on these three time window sets, which is represented as: Then iterate through For each time window in the data, for two adjacent time windows... , ,like If the time window is less than half an orbital period, then two adjacent time windows will be merged into one time window. ; Repeat the above process until... Within two adjacent time windows, there is no longer any In cases where the time window is less than half an orbital period, the final time window is obtained. This refers to the time window during which the optical satellite can see the target in the triangular area on the ground.

8. The method for calculating the visible window of an optical satellite for ground targets according to claim 7, characterized in that, Step S4 specifically includes: Let G be an arbitrary polygonal region on the ground. ,in The coordinates of the polygon vertices are the geodetic coordinates. Using the latitude and longitude of the vertices of the polygonal region as the horizontal and vertical coordinates respectively, the polygonal region is divided into sections using a planar polygon triangulation method. A triangular region, represented as in , , The coordinates of the vertices of the triangular region are given. By applying the aforementioned optical satellite visibility window algorithm for triangular regions on the ground, a set of time windows for each triangular region to be visible to the optical satellite can be obtained. This set is then combined to obtain the time window for the optical satellite to be visible to targets in any polygonal region on the ground, denoted as... 。

Images