Remote sensing satellite visible window calculation method based on geographic grid subdivision

Abstract

The invention discloses a geographic grid subdivision-based remote sensing satellite visible window calculation method, which comprises the following steps of: subdividing a global geographic grid to form a global discrete grid system; predicting to obtain a position vector sequence of the remote sensing satellite in a TEME reference system within a period of time based on an orbit prediction model, and converting the position vector sequence into an LLA coordinate system; mapping the sub-satellite point trajectory of the remote sensing satellite to a global discrete grid system to obtain a grid sequence of the sub-satellite point trajectory; mapping the geographic coordinates of the observation target to a global discrete grid system to obtain a target grid; based on the grid sequence and the target grid, obtaining a spherical distance sequence from the grid corresponding to the sub-satellite point to the target grid; and based on the spherical distance sequence and the observation load of the remote sensing satellite, calculating to obtain a visible window of the remote sensing satellite to the observation target. The method is applied to the field of space-based observation, and the calculation precision of the method can effectively meet the engineering requirements while the calculation efficiency of the visible window is high.

Description

Technical Field

[0001] The invention relates to the field of space-based observation technology, in particular to a remote sensing satellite visible window calculation method based on geographic grid division. Background Art

[0002] The visible window of a remote sensing satellite for ground targets refers to the time range in which the field of view of the payload it carries can observe the target. The calculation of the visible window is the basis and prerequisite for satellite earth observation mission planning. With the increase in the number of satellites in orbit and the development of satellite-related technologies, the rational allocation of satellite resources will become a new research trend. The problem of multi-satellite observation scheduling for multiple targets has gradually become a difficult and hot issue in research. Since the solution space of the scheduling problem is huge, and the calculation of the visible window is a necessary prerequisite for its solution, if the calculation efficiency of the visible window can be further improved, it will significantly shorten the calculation time of subsequent tasks such as observation mission planning, which has important research value and engineering significance.

[0003] The current visible window calculation method reads the ephemeris file line by line, calculates whether the geometric relationship between the satellite's spatial position and the target's geographical location at each moment satisfies the payload imaging constraint, and further selects the ephemeris moment corresponding to the appropriate roll angle and pitch angle as the imaging moment in the ephemeris that satisfies the constraint. Due to the lack of computing resources on board and the large number of ground targets, the calculation time is relatively long, which is suitable for high-performance equipment calculation on the ground, but difficult to meet the requirements of autonomous mission planning on board.

[0004] The earliest research on visibility issues focused on the visibility of satellites to point targets. The classic algorithm is the propagation algorithm, which determines visibility by traversing discrete time points in sequence with a fixed step length. This algorithm is accurate in calculation, but has low computational efficiency. It is often used as a comparative algorithm in the literature, and subsequent research goals are all focused on reducing computational time. For example, Lawton et al. first proposed using the visibility function method to analytically predict satellite-satellite and satellite-point target visibility, which greatly improved computational efficiency, but is only applicable to orbits with small eccentricity. Alfano et al. improved the adaptability of the visibility function so that it can be extended to any orbit; Ali et al. used a great circle arc to approximate the sub-satellite point trajectory under one orbital period. The model is simple and computationally efficient, but it is only applicable to low-orbit satellites; Mai et al. proposed a two-stage method of coarse search-fine estimation to quickly solve the visible window, which significantly shortened the calculation time; Wu et al. proposed a fast prediction method that can effectively shorten the imaging prediction time of satellite earth observation missions; Han et al. proposed an adaptive Hermite interpolation technology that can quickly determine the visibility of the satellite to the site, and the calculation time is greatly reduced; Zhang Jinxiu et al. proposed a fast prediction algorithm for the overpass and urban and rural area time of low-orbit satellites. The prediction accuracy of this algorithm can reach the second level, and it lays the foundation for the next step of precise search and prediction based on SGP4, which can greatly reduce the calculation amount of the on-board computer.

[0005] The above studies are all methods for calculating the visible window of satellites for point targets. For the visible resources of regional targets, the existing research methods can be divided into three categories: 1) Orbital period filtering: Most of the invisible periods are filtered out through satellite orbit and regional target information. For example, the document "Fast Algorithm for Time Window of Earth Observation Satellite Visiting Regional Targets" uses a great circle to approximate the sub-satellite point trajectory for conical and rectangular fields of view, which reduces the calculation time, but is only applicable to low-orbit satellites; 2) By reducing the discrete sampling points inside the region or at the region boundary, such as the reference "Fast Calculation Method of Time Window of Satellite to Regional Targets" considers the conical field of view and converts the problem into the union of the visible windows calculated by the time window of each boundary of the region; 3) Consider using parallel, two-stage search, dynamic step size and other computing methods to reduce the amount of calculation. For example, the literature "Real-time calculation method for satellite coverage time window" considers big data parallel calculation in the calculation process, and the literature "Dynamic step size fast algorithm for satellite time window calculation" and "A grid-based satellite access time window fast calculation method for regional targets" consider the dynamic step size in the calculation process. The literature "Fast calculation method for visibility of regional targets by remote sensing satellites" proposes to use binary search to increase the initial search step size, and only discretely sample the boundaries of regional targets to reduce the calculation scale. The above studies are all aimed at the calculation method of visible window of regional targets. The fastest time for simulating one day on ground computing equipment is about 0.1s, but there are still challenges in achieving efficient calculation on onboard computers.

[0006] The traditional visible time window calculation method determines visibility in sequence by traversing discrete time points with a fixed step length, and calculates the visible window of the payload based on the relative position relationship between the ground target and the satellite. This type of algorithm is accurate in calculation, but has high requirements for computing resources. At the same time, the transmission of ground calculation results to the satellite will increase the load of satellite-to-ground communication, causing delays in communication time, and it is difficult to meet the needs of emergency tasks. Satellite autonomous collaborative applications require on-board calculation of target visible windows to support satellite mission planning, and have high requirements for computing timeliness in emergency scenarios. At the same time, compared with the computing power of ground equipment, on-board computing resources are extremely limited, and often cannot meet the timeliness requirements in window calculations for emergency tasks or large-scale tasks. In addition, the satellite is more familiar with its own orbit and status, while the orbital information stored on the ground has errors, and the window accuracy of real-time calculation on the satellite is higher. Summary of the invention

[0007] In view of the high computational complexity and long computational time problems in the on-board autonomous visibility calculation of ground targets by conical and rectangular payload types of existing remote sensing satellites, the present invention provides a remote sensing satellite visibility window calculation method based on geographic grid division, which has high computational efficiency for the visibility window and its computational accuracy can effectively meet engineering requirements.

[0008] To achieve the above object, the present invention provides a method for calculating a remote sensing satellite visible window based on geographic grid subdivision, comprising the following steps: Step 1: divide the global geographic grid to form a global discrete grid system; Step 2, based on the orbit prediction model, obtain the position vector sequence of the remote sensing satellite in the TEM reference system within a period of time, and convert the position vector sequence into the LLA coordinate system; Step 3, based on the position vector sequence of the remote sensing satellite in the LLA coordinate system, mapping the sub-satellite point trajectory of the remote sensing satellite to the global discrete grid system to obtain a grid sequence of the sub-satellite point trajectory; Step 4, mapping the geographic coordinates of the observed target to the global discrete grid system to obtain the target grid; Step 5, based on the grid sequence and the target grid, obtaining a spherical distance sequence from the grid corresponding to the sub-satellite point to the target grid; Step 6: Based on the spherical distance sequence and the observation payload of the remote sensing satellite, the visible window of the remote sensing satellite for the observation target is calculated.

[0009] Compared with the prior art, the present invention has the following beneficial technical effects: Aiming at the characteristics of diversified remote sensing satellite observation tasks and uncertain geographical location, the present invention introduces a global discrete grid system to characterize the tasks, realizes standardized description of observation tasks, converts observation requirements into a set of predefined tasks, and realizes grid mapping of satellite resources based on a geographic grid partitioning framework and satellite orbit prediction technology. On the basis of grid characterization of observation tasks and satellite resources, the satellite-ground visible window is determined based on geometric relationship operations. Since the positional relationship between standardized grids can be pre-calculated and stored, the autonomous computing efficiency on board can be improved through data query. While having high computing efficiency for the visible window, its computing accuracy can effectively meet engineering requirements. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0011] Figure 1 It is a flow chart of a method for calculating a remote sensing satellite visible window based on geographic grid division in an embodiment of the present invention; Figure 2 It is a schematic diagram of a plane expansion of a regular icosahedron in an embodiment of the present invention; Figure 3 A schematic diagram of a hierarchical coding scheme in an embodiment of the present invention; Figure 4 Schematic diagram of the sixth level of the global discrete grid system in an embodiment of the present invention; Figure 5 Schematic diagram of the seventh level of the global discrete grid system in an embodiment of the present invention; Figure 6 Schematic diagram of visibility calculation of cone load mesh in an embodiment of the present invention; Figure 7 A schematic diagram of the visibility calculation of a rectangular load grid in an embodiment of the present invention; Figure 8 3 is an error diagram of a cone load with a semi-cone angle of 30° in an embodiment of the present invention.

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

[0013] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0014] It should be noted that all directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of the present invention are only used to explain the relative position relationship, movement status, etc. between the components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication will also change accordingly.

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

[0016] like Figure 1 The figure shows a method for calculating a visible window of a remote sensing satellite based on geographic grid division disclosed in this embodiment, which mainly includes the following steps: Step 1: divide the global geographic grid to form a global discrete grid system; Step 2: Based on the orbit prediction model, the position vector sequence of the remote sensing satellite in the TEM reference system within a period of time is obtained, and the position vector sequence is converted to the LLA coordinate system; Step 3, based on the position vector sequence of the remote sensing satellite in the LLA coordinate system, mapping the sub-satellite point trajectory of the remote sensing satellite to the global discrete grid system to obtain a grid sequence of the sub-satellite point trajectory; Step 4, mapping the geographic coordinates of the observed target to the global discrete grid system to obtain the target grid; Step 5, based on the grid sequence and the target grid, obtain the spherical distance sequence from the grid corresponding to the sub-satellite point to the target grid; Step 6: Based on the spherical distance sequence and the observation payload of the remote sensing satellite, the visible window of the remote sensing satellite for the observed target is calculated.

[0017] For global geographic grid division, the current mainstream methods include spherical triangle, spherical quadrilateral, spherical hexagonal grid and other division forms. Specifically: Spherical triangle grid is a spherical triangular mesh data structure that divides the earth's surface into a series of spherical triangles for analyzing spatial information. Triangulation is suitable for global terrain expression, has high data consistency and good data compression performance, and can reduce storage and transmission burdens. However, at the same time, the mesh division is relatively complex and the calculation is large, which is not suitable for mesh storage and fast encoding and decoding.

[0018] Spherical quadrilateral gridding is a method of dividing a sphere into quadrilateral regions. Its advantages are: 1) it has good geometric properties, which facilitates efficient spatial calculations; 2) it is suitable for data processing at different resolutions. However, quadrilaterals are prone to deformation at the poles, and the shape distortion is serious, which is not suitable for remote sensing satellite visibility calculations.

[0019] Spherical hexagonal gridding is a method of dividing the sphere into hexagonal regions. Its advantage is that the hexagonal regions have good distribution characteristics on the sphere, can better adapt to the changes on the earth's surface, and will not cause large deformation in the two-level regions, and improve the accuracy of data expression. At the same time, hexagons are more suitable for distance calculations between grids and for visible window calculations of remote sensing satellites for ground targets.

[0020] In this embodiment, a spherical hexagonal grid division method is adopted, and based on the division and related encoding methods of the open source system Uber-H3, a division and encoding scheme for the global discrete grid system is formed. In the specific implementation process, a projection method based on a regular polyhedron is adopted, and the regular icosahedron is used as the embedded regular polyhedron of the sphere, and the sphere is projected onto each face of the regular icosahedron. This is mainly because the icosahedron has the most faces compared to other Platonic cubes and fits the sphere best, so the projection distortion is minimal. The direction chosen for the regular icosahedron is Fuller's Dymaxion Orientation, which puts all twelve vertices of the regular icosahedron in the sea. The projection from the triangular plane of the regular icosahedron to the sphere is an inverse spherical center projection, and the tangent plane of the projection is parallel to the corresponding triangular plane of the regular icosahedron. Figure 2 This is the plane expansion diagram of a regular icosahedron. The grid formed by this division method is multi-level, with a total of 16 levels. The bottom layer (i.e., level 0) has 122 grids, including 12 regular pentagons, which are located at fixed points of the regular icosahedron, and 110 regular hexagons. Ignoring the pentagons at the vertices, this multi-level grid is subdivided using regular hexagons. The division aperture of the hexagonal grid is 7, that is, the area of ​​each parent grid is 7 times the area of ​​its child grid; the pentagonal grid has 5 hexagonal subgrids and 1 pentagonal subgrid. The advantage of hexagonal division is that the edge nearest neighbor, corner nearest neighbor, and distance nearest neighbor are consistent.

[0021] The encoding method of the global discrete grid system mainly adopts the hierarchical encoding method, providing a total of 15 levels of encoding hierarchy structure, and the index consists of 64 bits of binary. Figure 3 :The first bit, 6-8 bits are reserved and the value is 0; bits 2-5 represent the five index modes of the cell, bits 9-12 represent the corresponding coding level, ranging from [0,15]; bits 13-19 represent the 122 basic units of the global division when the coding level is 0, ranging from [0,121]; Figure 3In the , every three yellow digits represent the values ​​in the coordinate system of the face corresponding to coding level 1 to coding level 15, a total of 45 digits. When the coding level is 0, Uber-H3 divides the world into 122 basic units, including 110 hexagonal units and 12 pentagonal units centered on the vertices of the icosahedron, and assigns a digital index of [0,121] according to the latitude of the center point of the basic unit; basic unit 0 contains the northernmost center point, and basic unit 121 contains the southernmost center point. These units are selected to capture as many symmetries of the icosahedron as possible. Among them, the coordinate system uses IJK and FaceIJK coordinate systems for geographic grid coding. The IJK coordinate system refers to a coordinate system with three axes, X, Y, and Z, established on a hexagonal plane grid, and the interval between the three axes is 120°. Because the coding system is established on an icosahedron, in order to better encode, it is necessary to establish an IJK coordinate system on each face of the icosahedron, that is, to associate the face with the IJK coordinate system and establish the FaceIJK coordinate system.

[0022] Since most remote sensing satellites are deployed in low earth orbit (LEO), taking an orbit altitude of 500 kilometers as an example, the distance that its sub-satellite point trajectory moves on the earth's surface every second is about 6.7km, while the distance between the two grid centers of the sixth level of the global discrete grid system is about 5.6km. It can be seen that the sixth level grid can meet the resolution requirements of the sub-satellite point trajectory every second. The sub-satellite point trajectory of any LEO satellite is mapped to the sixth level grid every second, such as Figure 4 As shown. Since the distance between the two grid center points of the sixth level is very close to the distance that the LEO sub-satellite point trajectory moves on the earth's surface every second, this distance approximation makes the grids mapped by the sub-satellite point appear to be adjacently distributed for several consecutive seconds. From the perspective of spatial characteristics, the sixth-level grid has specific inherent properties such as size and direction; extending the consideration in the time dimension, when observing a longer period of time, based on the above grid characteristics, the grid sequence mapped by the sub-satellite point trajectory does not appear as a straight line, but is obviously composed of multiple short line segments. This phenomenon shows that the characteristics of the sixth-level grid have a significant impact on the grid mapping morphology of the sub-satellite point trajectory, making it difficult to form an ideal straight line distribution on a macro scale. The grids of the seventh level of the global discrete grid system are smaller than those of the sixth level, and the distance between adjacent grid centers is about 2.1km, which can better form an ideal straight line distribution, such as Figure 5 As shown. Its side length is about 1.22km, and for the mapping of point targets, it can also meet the resolution requirements in most scenarios. Therefore, in the specific implementation of steps 3 and 4, the sub-satellite point trajectory of the remote sensing satellite and the geographic coordinates of the observed target are mapped to the seventh-level grid in the global discrete grid system.

[0023] Given that the total number and spatial position of the seventh-level grids in the global discrete grid system are deterministic, the distance information between all grids can be pre-calculated and stored to facilitate real-time calculations on board. In practical applications, the onboard computer can directly read these stored data without repeating trigonometric functions and matrix operations for each observation task, thereby achieving rapid visible window determination. In specific applications, the geographic coordinates of the target and the satellite sub-satellite point are often expressed in the form of longitude and latitude. In the grid-based computing model, the longitude and latitude data need to be mapped to the corresponding grid code to achieve fast indexing, and the geographic grid division has an efficient encoding and decoding method. In the encoding process, the API function provided by H3 is used to convert the input geographic longitude and latitude into a hexagonal grid of a specified level containing the longitude and latitude point and the Hex2d coordinates corresponding to the grid. Decoding is the inverse process of encoding. The function receives the Hex2d coordinate encoding value corresponding to the input hexagonal grid and decodes the Hex2d coordinates into the geographic longitude and latitude data of the target point.

[0024] When the satellite is actually in operation, the satellite will autonomously maintain the satellite's ephemeris information in order to ensure mission completion and maintain attitude pointing. Compared with the calculation of ground equipment, the ephemeris data on the satellite will be more accurate, which can ensure higher-precision visible window calculation. The method in this embodiment first needs to implement the satellite orbit prediction and coordinate system conversion method, and on this basis, realize the visible window calculation of the conical load. It is worth noting that since the ephemeris will be maintained autonomously, no additional orbit prediction calculation is required when the visible window is determined autonomously on the satellite.

[0025] In this embodiment, SGP4 is used as the orbit prediction model, which takes into account factors such as the gravity of the earth, the gravity of the sun and the moon, atmospheric drag, and the perturbation caused by the rotation of the earth by analyzing the celestial mechanics equations, so as to achieve fast and accurate prediction of the satellite orbit. In specific implementation, the SGP4 model can obtain the initial orbital elements of the satellite, such as the semi-major axis, eccentricity, inclination, and ascending node longitude, by loading the classic TLE (Two-Line Element) data. The orbit prediction model SGP4 returns the original x, y, and z Cartesian coordinates of the satellite in the "True Equatorial Mean Equinox" (TEME) reference frame, which is an inertial reference frame (ECI) centered on the earth and does not rotate with the earth. However, SGP4 itself does not support the conversion of these coordinates to the more commonly used ECI frame (such as J2000 or ICRS), nor can it be directly converted to the Earth-centered Earth-fixed coordinate system (ECEF) or WGS84 coordinate system. The encoding and decoding of the grid requires the longitude and latitude information of the satellite, that is, the LLA coordinate system. To this end, it is necessary to perform coordinate transformation on the original output of the SGP4 model. Therefore, the specific implementation process of step 2 and step 3 in this embodiment is as follows: First, the TLE data of the remote sensing satellite is obtained and input into the orbit prediction module. The position vector sequence of the satellite in the TEM reference frame within a period of time is calculated through the SGP4 orbit prediction model. Subsequently, the position vector under the TEMES reference system is transformed into a coordinate system, and the position of the satellite is successively transformed from the TEMES reference system to the Earth-centered Earth-fixed (ECEF) reference system, and then further transformed into geographic coordinates (LLA) described by latitude, longitude and altitude; Furthermore, the obtained geographic coordinate sequence is mapped into a hexagonal grid and finally used as the input of the visibility calculation module; In the above process, orbit prediction is based on the satellite's maintained ephemeris data. Grid generation and distance calculation can be completed and stored in advance. Therefore, when performing autonomous visibility window calculation on board, the satellite only needs to complete the coordinate conversion, grid encoding and decoding, and visibility calculation of the orbit prediction results in real time. The coordinate conversion process is also required in the traditional visibility calculation method.

[0026] After completing the coordinate transformation of the satellite position vector sequence from the TEMES reference system to the LLA, the satellite's sub-satellite point trajectory can be further mapped to the hexagonal grid to generate a grid sequence P representing the satellite's sub-satellite point trajectory, which is:

[0027] in, t , λ t for t The latitude and longitude of the grid center corresponding to the subsatellite point at the moment, 1~T is the prediction period of the orbit prediction model, t= 1~T.

[0028] According to the longitude and latitude of the grid center mapped according to the observed target point 0. λ 0, the spherical distance sequence can be calculated , specifically:

[0029] in, d t for t The spherical distance from the grid corresponding to the sub-satellite point to the target grid at the moment, R For the radius of the earth, it can be taken as 6371km.

[0030] In the specific implementation process of step 6, when the observation load of the remote sensing satellite is a conical load, the calculation process of the visible window is: refer to Figure 6 , calculate the projection radius of the cone load on the earth's surfacer d ,for:

[0031] in, θ is the semi-cone angle of the conical load, h t for t The altitude of the remote sensing satellite at all times; The shortest distance After that, calculate the visibility sequence , specifically:

[0032] in, a t for t The visibility of the remote sensing satellite to the observation target at all times; The time corresponding to the subsequence of consecutive all-1s in the visibility sequence a is the visibility window.

[0033] In the specific implementation process of step 6, when the observation payload of the remote sensing satellite is a rectangular payload, the calculation process of the visible window is: refer to Figure 7 , calculate the projection distance of the maximum maneuvering range of the rectangular load in the side swing direction on the earth's surface q d ,for:

[0034] in, β is the half angle of the rectangular load in the roll direction, h t for t The altitude of the remote sensing satellite at all times; refer to Figure 7 , calculate the projection distance of the maximum maneuvering range of the rectangular load on the earth surface along the track direction p d ,for:

[0035] in, γ is the half angle of the rectangular load in the direction along the track; The shortest distance After that, calculate t Track distance at time b t ,for:

[0036] Compute visibility sequence , specifically:

[0037] in, a t for t The visibility of the remote sensing satellite to the observation target at all times; The time corresponding to the subsequence of consecutive all-1s in the visibility sequence a is the visibility window.

[0038] The following is a further explanation of the method for calculating the visible window of remote sensing satellites based on geographic grid division in this embodiment in conjunction with specific simulation experiments.

[0039] The start time of the test scenario is set to 04:00:00.000 on December 23, 2024, and the simulation time lasts for 1 hour, until 05:00:00.000 on December 23, 2024. The setting of this time period meets the verification requirements of the present invention for the visible window calculation of the satellite for the impending mission in the near time.

[0040] The design of the six-element parameters of the satellite is shown in Table 1. The satellite payload includes two types: conical payload and rectangular payload. The half-angle of the field of view of the payload can be adjusted according to user needs. In the simulation experiment, the half-angle of the conical payload was selected as 10°, 20°, 30°, 40° and 50°, and the vertical half-angle and horizontal half-angle of the rectangular payload field of view were set to 10°×10°, 20°×20°, 30°×30°, 40°×40° and 50°×50°. In order to verify the accuracy of the calculation results of the satellite's visibility to point targets, 15 targets were randomly selected near the subsatellite point trajectory to ensure that satellite payloads of different types of field of view sizes are visible to most targets, as shown in Table 2. Finally, the test results were compared and analyzed with the STK (Systems Tool Kit) software.

[0041] Table 1 Satellite orbit parameters

[0042] Table 2 Latitude and longitude of observation targets

[0043] Table 3 Calculation results of cone load with half cone angle of 30°

[0044] For 15 target points, when the semi-cone angle of the cone load is 30°, the comparison results with STK are shown in Table 3. From the results in Table 3, it can be seen that the average error between the start and end time of the window is less than 1s, which can meet the application requirements in engineering. The start and end time of the visible window of the cone load for 15 targets are compared with STK, and the error is less than 1.2s. Figure 8shown.

[0045] For the cone load, tests with different semi-cone angles are designed, and the results compared with STK are shown in Table 4. Since the geometric distance between the observation target and the sub-satellite point is realized by reading the pre-stored grid distance information, the visible window calculation based on geographic grid division can greatly reduce real-time mathematical operations, thereby improving efficiency and reducing calculation time. At the same time, it should be pointed out that due to the use of grids to approximate the positions of the observation target and the sub-satellite point, certain calculation errors will be caused. The average error between the start time and the end time of the window in the test results is less than 1s, and the calculation accuracy can meet the requirements of engineering applications.

[0046] Table 4 Comparison results of cone loads with different semi-cone angles

[0047] For the above 15 target points, when the field of view of the rectangular load is 30°×30°, the comparison results with STK are shown in Table 5. At the same time, for the rectangular load, experiments with different field sizes are designed, and the comparison results with STK are shown in Table 6. Among them, Target1 and Target2 are two randomly selected target points. It can be seen from the results in the table that the visible window calculation method based on geographic grid division, the window start time and end time are compared with the calculation results of STK, and the average error is in the order of 1s, indicating that the method proposed in this study can meet the application accuracy requirements of the project. It is also worth noting that the traditional visibility calculation method requires geometric relationship calculation for each target point. Therefore, as the number of tasks increases, the visible window calculation time will increase rapidly. When the number of targets is larger, the advantages of the method proposed in this embodiment will be more obvious. By storing the pre-calculated results, the complexity of real-time calculation can be greatly reduced.

[0048] Table 5 Calculation results of rectangular load field size 30°×30°

[0049] Table 6 Comparison results of rectangular loads with different field sizes

[0050] In summary, the maximum error between the visible window calculation method of remote sensing satellite based on geographic grid division in this embodiment and STK is about 1s in the visibility calculation of two different types of payloads, cone and rectangle, which can meet the requirements of engineering applications. It should also be pointed out that regional targets can also be characterized by grid division, and then the proposed method is used to calculate the visible window of the satellite to the target. After the grid division is completed, the calculation process is consistent with the point target calculation method.

[0051] The above description is only a preferred embodiment of the present invention, and does not limit the protection scope of the present invention. All equivalent structural changes made by using the contents of the present invention specification and drawings under the inventive concept of the present invention, or directly / indirectly applied in other related technical fields are included in the protection scope of the present invention.

Claims

1. A method for calculating the visible window of remote sensing satellites based on geographic grid division, characterized in that: The steps include: Step 1: divide the global geographic grid to form a global discrete grid system; Step 2, based on the orbit prediction model, obtain the position vector sequence of the remote sensing satellite in the TEM reference system within a period of time, and convert the position vector sequence into the LLA coordinate system; Step 3, based on the position vector sequence of the remote sensing satellite in the LLA coordinate system, mapping the sub-satellite point trajectory of the remote sensing satellite to the global discrete grid system to obtain a grid sequence of the sub-satellite point trajectory; Step 4, mapping the geographic coordinates of the observed target to the global discrete grid system to obtain the target grid; Step 5, based on the grid sequence and the target grid, obtaining a spherical distance sequence from the grid corresponding to the sub-satellite point to the target grid; Step 6: Based on the spherical distance sequence and the observation payload of the remote sensing satellite, the visible window of the remote sensing satellite for the observation target is calculated.

2. The method for calculating the visible window of remote sensing satellites based on geographic gridding according to claim 1, characterized in that: In step 1, a spherical hexagonal grid partitioning method is adopted, and based on the partitioning and related encoding methods of Uber-H3, a partitioning and encoding scheme of the global discrete grid system is formed.

3. The method for calculating the visible window of remote sensing satellites based on geographic grid division according to claim 2, characterized in that: In step 3 and step 4, the sub-satellite point trajectory of the remote sensing satellite and the geographic coordinates of the observed target are respectively mapped to the seventh-level grid in the global discrete grid system.

4. The method for calculating the visible window of remote sensing satellites based on geographic gridding according to claim 1, 2 or 3, characterized in that: In step 5, the spherical distance sequence , where 1~T is the prediction period of the orbit prediction model; Specifically: in, d t for t The spherical distance from the grid corresponding to the sub-satellite point to the target grid at the moment, R is the radius of the Earth, t , λ t for t The latitude and longitude of the grid center corresponding to the sub-satellite point at this moment, 0. λ 0 is the latitude and longitude of the target grid center.

5. The method for calculating the visible window of remote sensing satellites based on geographic grid division according to claim 4, characterized in that: In step 6, when the observation load of the remote sensing satellite is a conical load, the calculation process of the visible window is: Calculate the projected radius of the cone load on the earth's surface r d ,for: in, θ is the semi-cone angle of the conical load, h t for t The altitude of the remote sensing satellite at all times; The shortest distance After that, calculate the visibility sequence , specifically: in, a t for t The visibility of the remote sensing satellite to the observation target at all times; The time corresponding to the subsequence of consecutive all-1s in the visibility sequence a is the visibility window.

6. The method for calculating the visible window of remote sensing satellites based on geographic gridding according to claim 4, characterized in that: In step 6, when the observation payload of the remote sensing satellite is a rectangular payload, the calculation process of the visible window is: Calculate the projection distance of the maximum maneuvering range of the rectangular load in the side swing direction on the earth's surface q d ,for: in, β is the half angle of the rectangular load in the roll direction, h t for t The altitude of the remote sensing satellite at all times; Calculate the projection distance of the maximum maneuvering range of the rectangular load on the earth's surface along the track direction p d ,for: in, γ is the half angle of the rectangular load in the direction along the track; The shortest distance After that, calculate t Track distance at time b t ,for: Compute visibility sequence , specifically: in, a t for t The visibility of the remote sensing satellite to the observation target at all times; The time corresponding to the subsequence of consecutive all-1s in the visibility sequence a is the visibility window.

Images