A satellite access information service method based on global grid

Abstract

The present invention discloses a satellite access information service method based on global grid division, which belongs to the field of satellite remote sensing. Based on the theory of global grid division, the present invention divides the world into several grids, which seamlessly cover the world without intersection, calculate the satellite's access information to each grid in advance and store it; when a remote sensing observation demand is received, the corresponding grid of the observation demand is quickly calculated according to the observation demand position, and the satellite's access information to the grid is retrieved from the database, and it is used as the satellite's access information to the observation demand. The present invention adopts a method of exchanging space for time, and has the characteristics of high computing efficiency, easy implementation, easy expansion, and strong compatibility. It can solve the problem that large-scale remote sensing constellations take a long time to access and calculate massive ground observation demands and cannot meet timeliness requirements.

Description

Technical Field

[0001] The invention belongs to the field of satellite remote sensing, and in particular refers to a satellite access information service method based on a global grid, which can be used for rapid calculation of access to massive observation requirements by large-scale remote sensing star clusters in the field of satellite remote sensing mission planning. Background Art

[0002] In the field of remote sensing satellite mission planning, calculating the satellite's access information such as the access time and observation angle for observation needs is the basis for all subsequent work such as satellite mission planning and instruction compilation. The traditional method is the real-time calculation mode, which calculates the access information for observation needs point by point according to the ephemeris file. Since the ephemeris file is large, the calculation time is long and the calculation efficiency is low.

[0003] In order to improve computing efficiency, the project uses multi-node parallel computing to improve the overall efficiency of multiple tasks. On the other hand, it optimizes the computing method to reduce invalid computing of a single node and further improve computing efficiency. However, these methods are still real-time computing modes. In the context of limited growth in the amount of tasks and a controllable number of satellites, this mode can ensure timeliness requirements by expanding server resources. However, when faced with the massive observation needs of large-scale remote sensing constellations, the amount of tasks grows exponentially, and server resources are difficult to expand infinitely, making it difficult to ensure high timeliness requirements. In order to meet the high timeliness computing requirements of massive tasks, it is urgent to study a new technical system to achieve a leapfrog improvement in efficiency. Summary of the invention

[0004] In view of this, the present invention proposes a satellite access information service method based on global grid. Aiming at the problem that the access calculation of massive ground observation needs of future large-scale remote sensing constellations takes a long time and cannot meet the timeliness requirements, the present invention changes the real-time calculation into data query, which can greatly improve the operation efficiency and has great application value.

[0005] The technical solution adopted by the present invention is:

[0006] A satellite access information service method based on global gridding comprises the following steps:

[0007] (1) Divide the earth's surface into seamless, non-overlapping, and scale-continuous multi-level grids based on longitude and latitude coordinates, store grid access information, and each grid has a unique code. All grids have similar shapes, and the lower grids are obtained by splitting the upper grids.

[0008] (2) Divide the ephemeris into multiple segments, and determine whether the period meets the lighting conditions required for optical remote sensing satellite observation. The period that does not meet the lighting conditions is directly removed and no longer participates in the calculation;

[0009] (3) Set the time step and the side view step. For the time period that meets the lighting conditions, calculate the access information of each grid that the satellite can cover in each time period according to the two dimensions of time step and side view step. The access information includes the access time and the side view angle. If the coverage area of ​​a time period partially meets the imaging conditions, calculate whether each grid meets the imaging conditions and save the access information of the grids that meet the imaging conditions into the database.

[0010] (4) Calculate the code of the grid where the observation requirement is located based on the longitude and latitude of the observation requirement, and retrieve the satellite's access information to the grid from the database as the satellite's access information to the observation requirement.

[0011] Furthermore, the specific method of step (1) is:

[0012] (101) The grid level n is determined according to the ground width observed by satellite, where n ≥ 2. The world is divided into a multi-level grid system according to longitude and latitude. The first level division method is:

[0013] Horizontally, from the equator to the north and south poles, each latitude difference of 6° is a row, 15 rows in each hemisphere, a total of 30 rows; vertically, from longitude 0° to the east and west, each longitude difference of 6° is a column, 30 columns in each hemisphere, a total of 60 columns;

[0014] Divide the first-level grid into 4 equal parts to obtain the second-level grid, and divide the second-level grid into 4 equal parts to obtain the third-level grid, and so on, until the n-th level grid;

[0015] (102) Construct a grid coding rule, including layer number, northern and southern hemisphere code, row number, eastern and western hemisphere number, and column number, where:

[0016] The layer number is a positive integer, the first layer is 1, the second layer is 2, and so on;

[0017] The codes for the northern and southern hemispheres are: N for the northern hemisphere, S for the southern hemisphere;

[0018] The row number is a positive integer, starting from 1, and increasing by 1 from the equator to the second level. The maximum value of the row number is different for different levels. The maximum value of the first level is 15, and it doubles with each level.

[0019] The Eastern and Western Hemispheres are numbered as follows: E for the Eastern Hemisphere and W for the Western Hemisphere;

[0020] The column number is a positive integer, starting from 1, and increasing by 1 from the 0 degree meridian to both sides. The maximum value of the column number is different at different levels. The maximum value of the first level is 30, and it doubles at each level.

[0021] (103) Calculate the longitude and latitude range of the grid according to the grid number; if the layer number of the grid is lid, the row number is row, and the column number is col, then the longitude and latitude range of the grid is calculated as follows:

[0022] The latitude and longitude difference of this layer grid is Δ = 6 / 2 lid-1 ;

[0023] Grid minimum longitude minLon = (col - 1) × Δ;

[0024] Grid maximum longitude maxLon = col×Δ;

[0025] Grid minimum latitude minLat = (row - 1) × Δ;

[0026] Grid maximum latitude maxLat = row×Δ;

[0027] If the code for the eastern and western hemispheres is W, let minLon = -minLon, maxLon = -maxLon;

[0028] If the code for the northern and southern hemispheres is S, let minLat = -minLat, maxLat = -maxLat;

[0029] (104) According to the longitude and latitude coordinates (lon, Lat) and the layer number lid, calculate the grid number where it is located:

[0030] The latitude and longitude difference of this layer grid is Δ = 6 / 2 lid-1 ;

[0031] Grid column number col = TopInt(lon / Δ)

[0032] Grid row number row = TopInt(Lat / Δ)

[0033] TopInt is the upward rounding function;

[0034] If lon ≥ 0, the code for the eastern and western hemispheres is E, otherwise, the code for the eastern and western hemispheres is W;

[0035] If Lat ≥ 0, the code for the northern and southern hemispheres is N, otherwise, the code for the northern and southern hemispheres is S;

[0036] (105) A grid layer with a grid width at the equator that matches the satellite observation width is used as the calculation carrier. The first layer of grids has 60 grids per row. The width of the first layer of grids at the equator is 666 km. The width of each subsequent layer of grids at the equator is half of the width of the previous layer of grids at the equator and is rounded up.

[0037] (106) A two-layer storage architecture is used to store grid access information. The first layer contains a table named SatGridLevel, which stores the applicable grid levels of different satellites and is used to quickly determine the grid level corresponding to the satellite during grid access calculation;

[0038] The second layer contains multiple tables, which store the satellite's access information to the grid. Each satellite corresponds to a table, and the table name is: satellite code + AI.

[0039] Furthermore, the specific method of step (2) is:

[0040] (201) The time period of the segmented ephemeris is set to 2 minutes, and each segment of the ephemeris is a set of subsatellite points SEG = {p1, p2, ..., p i-1 ,p i ,p i+1 ,……,p 120}, for any subsatellite point p in SEG i , if any of the following conditions exist:

[0041]

[0042]

[0043] Then p i The points after that are separated from SEG, generating a new segmented ephemeris SEG1={p i ,p i+1 ,……,p 120};

[0044] (202) Using the satellite pointing calculation method, calculate each segment of the ephemeris {p1, p2, ..., p n}, the specific method is:

[0045] The earth surface equation in the J2000 coordinate system is:

[0046] (Formula 1)

[0047] Among them, a and b are the equatorial radius of the earth, a=b=6378137m;

[0048] c is the polar radius of the Earth, c=a×(1-f);

[0049] f is the earth's flattening, f=1 / 298.257223563;

[0050] Satellite boresight vector in J2000 coordinate system The parametric equation of the straight line vector is:

[0051] (Formula 2)

[0052] Among them, x0, y0, and z0 are the spatial positions of the satellite in the J2000 coordinate system at a certain time, which are obtained from the ephemeris file; when the satellite side viewing angle is α, m, n, and p are the satellite line of sight vectors In the three components of the coordinate axis, the specific calculation process is as follows:

[0053] Satellite boresight in orbital coordinate system ; (Formula 3)

[0054] Satellite boresight in J2000 coordinate system =(m,n,p)= ; (Formula 4)

[0055] in, The transformation matrix from the orbital system to the J2000 coordinate system. The elements of this transformation matrix are:

[0056] (Formula 5)

[0057] Among them, i=1,2,3 corresponds to the three components of each row vector in the transformation matrix C;

[0058] C(1,i), C(2,j) and C(3,k) are the row vectors of the first, second and third rows of the matrix C respectively;

[0059] Represents the position vector of the satellite in the J2000 coordinate system, obtained from the ephemeris file;

[0060] Represents the velocity vector in the satellite J2000 coordinate system, obtained from the ephemeris file;

[0061] Combine (Formula 1)-(Formula 4) and solve the equations to calculate the intersection of the satellite's visual axis and the earth's surface, that is, the ground pointing point of the satellite at a certain time and with a certain field of view angle;

[0062] In each ephemeris segment, the satellite's coverage of the ground is a quadrilateral, and the vertices of the quadrilateral are (P1, P2, P3, P4), where:

[0063] Substitute the satellite J2000 spatial position and the minimum side viewing angle range at the beginning of each ephemeris into (Formula 1)-(Formula 4) to calculate the P1 point;

[0064] Substitute the satellite J2000 spatial position and the maximum side viewing angle range at the beginning of each ephemeris into (Formula 1)-(Formula 4) to calculate the P2 point;

[0065] Substitute the satellite J2000 spatial position and minimum side viewing angle range at the last moment of each ephemeris into (Formula 1)-(Formula 4) to calculate the P3 point;

[0066] Substitute the satellite J2000 spatial position and the maximum side viewing angle range at the last moment of each ephemeris into (Formula 1)-(Formula 4) to calculate the P4 point;

[0067] (203) If the satellite observation has requirements for lighting conditions, then determine the lighting conditions of the segmented ephemeris coverage range, and calculate the solar altitude angles sun1, sun2, sun3, and sun4 of points P1, P2, P3, and P4 respectively. If the four solar altitude angles are all less than 0, the area does not meet the observation conditions and is directly discarded; if the four solar altitude angles are all greater than or equal to 0, the area fully meets the imaging conditions; if there are K solar altitude angles greater than 0, K = 1, 2, 3, then the area partially meets the imaging conditions;

[0068] The calculation formula for the solar altitude angle is:

[0069] (Formula 6)

[0070] Where: θ — the solar altitude angle at the observation point, in degrees;

[0071] — latitude of the observation site, in degrees;

[0072] δ — solar declination, in degrees;

[0073] τ — observation geotime angle, in degrees;

[0074] The calculation formula for the solar declination δ is:

[0075] δ = 23.45 × sin [(360 / 365) × (284 + n)] (Formula 7)

[0076] Where n is the number of days in a year, calculated from the time in the ephemeris;

[0077] The calculation formula of the observation time angle τ is:

[0078] (Formula 8)

[0079] in:

[0080] st is the true solar time, that is, the local time required for observation. For every 15 degrees of longitude on the earth's surface, the local time differs by 1 hour;

[0081] Δt represents the time difference between the local time of the observation requirement and the Coordinated Universal Time;

[0082] lon is the longitude of the observation requirement;

[0083] t is the corresponding ephemeris time;

[0084] tz is the time zone number corresponding to the ephemeris time standard. If the ephemeris time standard is Coordinated Universal Time, tz=0; if the ephemeris time standard is Beijing Time, tz=8;

[0085] (204) Save the access information of the time period into the segmented filtering information list: SegList={s i |i=1,2,…}, where s i Contains the start time t1, end time t2 of the time period, and whether the imaging conditions are met.

[0086] Furthermore, the specific method of step (3) is:

[0087] (301) Select the satellite field of view angle as the step length of the side viewing angle to ensure that the satellite pointing crosses one grid each time the side viewing angle steps. When the satellite is in a high latitude area, the distance in the longitude direction becomes smaller, the width of the grid becomes smaller, and the side viewing angle steps once. If the satellite pointing crosses multiple grids and these grids are all within the range of the satellite field of view angle Vα, then the side viewing angles of these grids are the same.

[0088] (302) Set the time step length to the ratio of the satellite observation width Vw to the distance between adjacent sub-satellite points, that is, the time step length ts = Vw / 7.8; if the satellite observation width Vw is less than 7.8, set ts = 1. In this case, if the satellite crosses multiple grids during one time step, the access time of these grids is still considered to be the same;

[0089] (303) Traverse the pre-screening time period set SegList that meets the conditions, for any time period s i 、 Satellite side viewing range (α min , α max ), and process as follows:

[0090] ①Let t = t1, α = α min ;

[0091] ② When t>t2, execute step (304);

[0092] ③When α>α max , execute step ⑧;

[0093] ④ Calculate the ground pointing point p of the satellite at time t and side angle α;

[0094] ⑤If s iIf the imaging conditions of the time period are fully met, then execute step ⑥; if they are partially met, then calculate the solar altitude angle of point p at time t; if the solar altitude angle does not meet the imaging conditions, then execute step ⑦;

[0095] ⑥ Calculate the grid number Gid where point p is located according to the method of step (104), record the grid access information GAI, set GAI.gid=Gid, GAI.t=t, GAI.α=α, and save the grid access information to the grid access information list GAIList;

[0096] ⑦Let α=α+Vα, and execute step ③;

[0097] ⑧Let t=t+ts, and execute step ②;

[0098] (304) For high-latitude areas and grid widths less than 7.8 km, the grid access information of GAIList is checked. When the grid column number in a row is crossed, the access information of the crossed grid is supplemented, and the access information of the crossed grid is equal to the access information of the grid in the previous column in the row; when the grid row number in a column is crossed, the access information of the crossed grid is supplemented, and the access information of the crossed grid is equal to the access information of the grid in the previous row in the column;

[0099] (305) The access information is saved in the database; in order to improve the efficiency of data retrieval, outdated data is deleted, that is, data with an access time earlier than the current time is deleted.

[0100] Furthermore, the specific method of step (4) is:

[0101] (401) After receiving the observation requirement, according to the resolution requirement in the observation requirement, obtain the applicable satellite, determine which grid access information table to query next, and put it into the table set to be queried TabList;

[0102] (402) According to the observation requirement location, the grid number Gid corresponding to the observation requirement is calculated in the manner of step (104);

[0103] (403) According to the grid number Gid, each table in TabList is queried to obtain the access satellite, access time and access angle required for observation.

[0104] Compared with the calculation method in the background art, the present invention has the following advantages:

[0105] 1. The core idea of ​​the present invention is to trade space for time, and change real-time calculation into data query, which can greatly improve the operation efficiency and has great application value. Specifically, based on the theory of global grid division, the present invention divides the world into a multi-level grid system, with different grid sizes at different levels, and each grid has seamless coverage of the world without intersection; the applicable grid level of each satellite is determined according to the satellite observation parameters, and the satellite's access information to the applicable grid level is calculated in advance and stored; when a remote sensing observation demand is received, the grid where the observation demand is located is quickly calculated according to the observation demand location, and the satellite's access information to the grid is retrieved from the database, and it is used as the satellite's access information to the observation demand.

[0106] 2. Compared with the conventional algorithms in the background technology, the present invention proposes a new working mechanism, which changes the original real-time calculation mode into an advance calculation and real-time query mode, and exchanges the storage space of the advance calculation results for the time consumed by the real-time calculation. Although the time consumption of the advance calculation is comparable to that of the conventional method, it greatly improves the real-time response capability of the entire system tasks.

[0107] 3. The parallel computing architecture in the background technology and the idea of ​​reducing the amount of computing through optimization methods are still applicable. The present invention has a good inheritance relationship with the original computing architecture and computing resources, but only expands the storage service architecture of the computing results based on the original computing architecture.

[0108] 4. The present invention has the characteristics of simple software architecture, easy implementation and debugging, easy expansion and strong compatibility. BRIEF DESCRIPTION OF THE DRAWINGS

[0109] Figure 1 This is a schematic diagram of the global grid division method.

[0110] Figure 2 It is a flow chart of segmented screening process.

[0111] Figure 3 It is a flow chart of grid access information calculation.

[0112] Figure 4 It is a schematic diagram of the grids that are crossed during grid access calculation.

[0113] Figure 5 It is the flow chart of observation demand access information service. DETAILED DESCRIPTION

[0114] Below, the present invention is further described in conjunction with the accompanying drawings.

[0115] A satellite access information service method based on global gridding comprises the following steps:

[0116] (1) Global grid architecture design: Divide the earth's surface into multi-level grids that are seamless, non-overlapping, and continuous in scale based on longitude and latitude coordinates. All grids are similar in shape, and the lower grids can be obtained by splitting the upper grids. Each grid has a unique code. The grid level and grid latitude and longitude range can be calculated based on the grid code. The grid code can be quickly calculated based on the geographic coordinates and the specified grid level. Determine the applicable grid level based on the ground width observed by satellites. Design a storage architecture for grid access information.

[0117] (2) Time period pre-screening: A large number of optical remote sensing satellite observations require good lighting conditions. Before calculating grid access information, the ephemeris is divided into several segments, and the time period that does not meet the lighting conditions is judged segment by segment. The time period that does not meet the lighting conditions is directly removed and no longer participates in the calculation, which can greatly reduce the amount of calculation and improve the overall calculation efficiency.

[0118] (3) Calculation of grid access information: For the time period that meets the imaging conditions, the satellite access information for each grid that can be covered during the time period is calculated according to a certain time step and side viewing angle step. The information mainly includes the access time and side viewing angle. If the area covered during the time period partially meets the imaging conditions, it is necessary to calculate whether each grid meets the imaging conditions and save the access information of the grids that meet the imaging conditions into the database.

[0119] (4) Observation demand access information service: The code of the grid where the observation demand is located is calculated based on the longitude and latitude of the observation demand, and the satellite access information of the grid is retrieved from the database and used as the satellite access information for the observation demand.

[0120] The specific implementation of step (1) is as follows:

[0121] (101) Grid division method: The world is divided into a multi-level grid system based on longitude and latitude. The first-level division method is as follows: horizontally, from the equator to the north and south poles, each 6° difference in latitude is a row, 15 rows in each hemisphere, a total of 30 rows; vertically, from longitude 0° to the east and west, each 6° difference in longitude is a column, 30 columns in each hemisphere, a total of 60 columns; the first-level grid is divided into 4 equal parts to obtain the second-level grid, and the second-level grid is further divided into 4 equal parts to obtain the third-level grid, and so on, it can be infinitely subdivided, such as Figure 1 shown.

[0122] (102) Grid coding rules: layer number + northern and southern hemisphere code + row number + eastern and western hemisphere number + column number. Among them:

[0123] The layer number is a positive integer, the first layer is 1, the second layer is 2, and so on;

[0124] The codes for the northern and southern hemispheres are: N for the northern hemisphere, S for the southern hemisphere;

[0125] The row number is a positive integer, starting from 1, and increasing by 1 from the equator to the second level. The maximum value of the row number is different for different levels. The maximum value of the first level is 15, and it doubles with each level.

[0126] The Eastern and Western Hemispheres are numbered as follows: E for the Eastern Hemisphere and W for the Western Hemisphere;

[0127] The column number is a positive integer, starting from 1, and increasing by 1 from the 0 degree meridian to both sides. The maximum value of the column number is different at different levels. The maximum value of the first level is 30, and it doubles at each level.

[0128] (103) Calculate the longitude and latitude range of the grid according to the grid number: The layer number of the grid is lid, the row number is row, and the column number is col. The longitude and latitude range of the grid is calculated as follows:

[0129] The longitude and latitude difference of this layer grid Δ =6 / 2 lid-1 ;

[0130] Grid minimum longitude minLon = (col - 1) × Δ;

[0131] Grid maximum longitude maxLon = col×Δ;

[0132] Grid minimum latitude minLat = (row - 1) × Δ;

[0133] Grid maximum latitude maxLat = row×Δ;

[0134] If the code for the eastern and western hemispheres is W, then minLon = -minLon, maxLon = -maxLon;

[0135] If the code for the northern and southern hemispheres is S, then minLat = -minLat, maxLat = -maxLat.

[0136] (104) According to the longitude and latitude coordinates (lon, Lat) and the layer number lid, calculate the grid number where it is located:

[0137] The longitude and latitude difference of this layer grid Δ =6 / 2 lid-1 ;

[0138] Grid column number col = TopInt(lon / Δ)

[0139] Grid row number row = TopInt(Lat / Δ)

[0140] TopInt is the upward rounding function;

[0141] If lon ≥ 0, the code for the eastern and western hemispheres is E, otherwise, the code for the eastern and western hemispheres is W;

[0142] If Lat ≥ 0, the northern and southern hemispheres are coded as N, otherwise, the northern and southern hemispheres are coded as S.

[0143] (105) Determine the applicable grid level: In order to ensure the accuracy of satellite calculation of observation demand access information, a grid level with a grid width at the equatorial position that matches the satellite observation width is used as the calculation carrier. The first-level grid has 60 grids per row. The earth's equator is about 40,000 km long. The width of the first-level grid at the equator is 40,000 / 40≈666 km. The second-level grid is 4 equal parts of the first-level grid. The width of the second-level grid at the equator is 666 / 2=333 km, and so on. The width of each level of grid at the equator is shown in Table 1 below:

[0144]

[0145] (106) Grid access information storage architecture design: Remote sensing satellites have various resolutions. The ground observation width of high-resolution satellites is relatively small, usually below 10 km. Satellites can access a large number of grids every day. In order to improve service efficiency, a two-layer storage architecture is used for grid access information. The first layer contains a table named SatGridLevel, which stores the applicable grid levels of different satellites. It is used to quickly determine the grid level corresponding to the satellite during grid access calculation. The table structure is shown in Table 2 below:

[0146]

[0147] The second layer contains multiple tables to store satellite access information to the grid. Each satellite has a table. The table name is: satellite code + AI. The table structure is shown in Table 3 below:

[0148] .

[0149] The processing flow of step (2) is as follows Figure 2 As shown, the specific method is:

[0150] (201) Taking into account the east-west span of the segmented ephemeris and the amount of pre-screening calculation, the time period of the segmented ephemeris is set to 2 minutes. Each segment of the ephemeris is a set of subsatellite points SEG = {p1, p2, ..., p i-1 ,p i ,p i+1 ,……,p 120}, for the convenience of calculation, it is necessary to ensure the continuity of the longitude and latitude coordinates of the sub-satellite point. The specific calculation method is as follows:

[0151] For any subsatellite point p in SEG i,If any of the following conditions exist, then pi and subsequent points are separated from SEG to generate a new segmented ephemeris SEG1={p i ,p i+1 ,……,p 120}:

[0152]

[0153]

[0154] (202) Using the satellite pointing calculation method, calculate each segment of the ephemeris {p1, p2, ..., p n}Satellite coverage, satellite pointing calculation method is:

[0155] The earth surface equation in the J2000 coordinate system is:

[0156] (Formula 1)

[0157] Where a and b are the equatorial radii of the earth, a=b=6378137m;

[0158] c is the polar radius of the Earth, c=a×(1-f);

[0159] f is the earth's flattening, f=1 / 298.257223563;

[0160] Satellite boresight vector in J2000 coordinate system The parametric equation of the straight line vector is:

[0161] (Formula 2)

[0162] x0, y0, z0 are the spatial positions of the satellite in the J2000 coordinate system at a certain time, which can be obtained from the ephemeris file; when the satellite side view angle is α, m, n, p are the satellite line of sight vectors In the three components of the coordinate axis, the specific calculation process is as follows:

[0163] Satellite line of sight in orbital coordinate system: ; (Formula 3)

[0164] Satellite boresight in J2000 coordinate system: =(m,n,p)= ; (Formula 4)

[0165] The transformation matrix from the orbital system to the J2000 coordinate system. The elements of this transformation matrix are:

[0166] (Formula 5)

[0167] Among them: i=1,2,3 corresponds to the three components of each row vector in the transformation matrix C;

[0168] C(1,i), C(2,j) and C(3,k) are the row vectors of the first, second and third rows of the matrix C respectively;

[0169] Represents the position vector of the satellite in the J2000 coordinate system, which can be obtained from the ephemeris file;

[0170] Represents the velocity vector in the satellite J2000 coordinate system, which can be obtained from the ephemeris file.

[0171] By solving equations (1)-(4) simultaneously, we can calculate the intersection of the satellite's line of sight and the earth's surface, that is, the point on the ground where the satellite is pointing at a certain field of view at a certain moment.

[0172] In each ephemeris segment, the satellite's coverage of the ground is a quadrilateral, and the vertices of the quadrilateral are (P1, P2, P3, P4), where:

[0173] The calculation process of P1 point is to substitute the satellite J2000 spatial position and the minimum side viewing angle range at the beginning of each ephemeris into (Formula 1)-(Formula 4) to calculate;

[0174] The calculation process of P2 point is to substitute the satellite J2000 spatial position and the maximum side viewing angle range at the beginning of each ephemeris into (Formula 1)-(Formula 4) to calculate;

[0175] The calculation process of P3 point is to substitute the satellite J2000 spatial position and the minimum side viewing angle range at the last moment of each ephemeris into (Formula 1)-(Formula 4) to calculate;

[0176] The calculation process of P4 point is to substitute the satellite J2000 spatial position and the maximum side viewing angle range at the last moment of each ephemeris into (Formula 1)-(Formula 4) to calculate it.

[0177] (203) If the satellite observation has requirements for lighting conditions, then determine the lighting conditions of the segmented ephemeris coverage area and calculate the solar altitude angles (sun1, sun2, sun3, sun4) of (P1, P2, P3, P4) respectively. If the four solar altitude angles are all less than 0, the area does not meet the observation conditions and is directly discarded; if the four solar altitude angles are all greater than or equal to 0, the area fully meets the imaging conditions; if 1-3 solar altitude angles are greater than 0, the area partially meets the imaging conditions.

[0178] The formula for calculating the solar altitude angle is:

[0179] (Formula 6)

[0180] Where: θ — the solar altitude angle at the observation point, in degrees;

[0181] — latitude of the observation site, in degrees;

[0182] δ — solar declination, in degrees;

[0183] τ — Observation geotime angle, in degrees.

[0184] The approximate calculation formula for the solar declination δ is:

[0185] δ = 23.45 × sin [(360 / 365) × (284 + n)] (Formula 7)

[0186] Where n is the number of days in a year, which can be calculated from the time in the ephemeris.

[0187] The calculation formula of the observation time angle τ is:

[0188] (Formula 8)

[0189] Wherein, st is the true solar time, i.e. the local time required for observation. The local time differs by 1 hour for every 15 degrees of longitude on the earth's surface; Δt represents the time difference between the local time required for observation and the coordinated universal time; lon is the longitude required for observation; t is the corresponding ephemeris time; tz is the time zone number corresponding to the ephemeris time standard. If the ephemeris time standard is the coordinated universal time, tz=0; if the ephemeris time standard is Beijing time, tz=8;

[0190] (204) The access information of this period is saved in the segmented filtering information list: SegList={s i |i=1,2,…}, s i The main information is the start time, end time, and imaging condition satisfaction of the period (i.e., "fully satisfied" or "partially satisfied").

[0191] The processing flow of step (3) is as follows Figure 3 As shown, the specific method is:

[0192] (301) Determine the side viewing angle step length: Because the grid level is determined by the satellite observation width and the grid width near the equator, the satellite observation width is determined by the satellite altitude and the field of view angle Vα. Therefore, the satellite field of view angle is selected as the side viewing angle step length to ensure that the satellite pointing will cross one grid each time the side viewing angle steps forward. When the satellite is in a high-latitude area, the distance in the longitude direction becomes smaller, the grid width becomes smaller, and the satellite pointing may cross multiple grids each time the side viewing angle steps forward, but these grids are all within the range of the satellite field of view angle Vα, that is, the side viewing angles of these grids are the same.

[0193] (302) Determine the time step length: The satellite flying speed is about 7.8 km / s, and the distance between adjacent sub-satellite points with a 1-second interval ephemeris is also about 7.8 km. To ensure that the satellite pointing crosses one grid in one time step, the time step length should be the ratio of the equatorial grid width (satellite observation width Vw) to the distance between adjacent sub-satellite points, that is, the time step length ts = Vw / 7.8. If the satellite width Vw is less than 7.8, let ts = 1. In this case, the satellite may cross multiple grids in one time step. Since the satellite observation time is greater than 1 second, it can be considered that the access time of these grids is the same.

[0194] (303) Traverse the pre-screening time period set SegList that meets the conditions, for any time period s i (t1, t2, imaging conditions), satellite side viewing range (α min , α max ), and process as follows:

[0195] ①Let t = t1, α = α min ;

[0196] ② When t>t2, execute step (304);

[0197] ③When α>α max , execute step ⑧;

[0198] ④ Calculate the ground pointing point of the satellite at time t with an angle of view from the α side;

[0199] ⑤If s i If the imaging conditions of the time period are fully met, execute step ⑥. If they are partially met, calculate the solar altitude angle of point p at time t. If the solar altitude angle does not meet the imaging conditions, execute step ⑦.

[0200] ⑥ According to step (104), the grid number Gid where point p is located is calculated, and the grid access information GAI is recorded, and GAI.gid=Gid, GAI.t=t, GAI.α=α, and the grid access information is saved in the grid access information list GAIList;

[0201] ⑦ If α=α+Vα, execute step ③;

[0202] ⑧If t=t+ts, execute step ②;

[0203] (304) Grid access information check and interpolation, mainly for high latitudes and grid width less than 7.8 km. In high latitudes, a side view step may span multiple columns of grids. When the grid width is less than 7.8 km, a time step may span multiple rows of grids, such as Figure 4 shown.

[0204] Check the grid access information of GAIList. When the grid column number in a row is crossed, the access information of the crossed grid needs to be supplemented. The access information of the crossed grid is equal to the access information of the grid in the previous column in this row. When the grid row number in a column is crossed, the access information of the crossed grid needs to be supplemented. The access information of the crossed grid is equal to the access information of the grid in the previous row in this column.

[0205] (305) Save the past access information to the database, and the table structure is shown in Table 3. In order to improve the efficiency of data retrieval, it is necessary to delete the outdated data, that is, delete the data whose access time is earlier than the current time.

[0206] The processing flow of step (4) is as follows Figure 5 As shown, the specific method is:

[0207] (401) After receiving the observation request, according to the resolution requirement in the observation request, query the satellite applicable grid level table (as shown in Table 2), obtain the applicable satellite, the corresponding grid level, and the grid access information table (as shown in Table 3) set, and put it into the to-be-queried table set TabList;

[0208] (402) according to the observation requirement location and the applicable grid level queried in step (401), the grid number Gid corresponding to the observation requirement is calculated according to step (104);

[0209] (403) According to the grid number Gid, each table in TabList is queried to obtain all access information (satellite code, access time, side view angle) of the grid. All access information is the access information of the requirement. There may be multiple access information corresponding to the observation requirement.

[0210] In summary, the present invention is based on the global grid theory, divides the world into a multi-level grid system, calculates the access information of each satellite to the grid in advance, and stores it in the database; after receiving the observation demand, calculates the grid number according to its position, and then queries the access information of the grid from the database as the access information of the demand, converts the original real-time calculation mode into a data query service, greatly improves the system operation efficiency, and is particularly suitable for the access information calculation scenario of large-scale star clusters for massive observation needs. The present invention has the characteristics of high operation efficiency, easy implementation, easy expansion, and strong compatibility.

Claims

1. A satellite access information service method based on global gridding, characterized in that: The following steps are involved: (1) Divide the earth's surface into seamless, non-overlapping, and scale-continuous multi-level grids based on longitude and latitude coordinates, store grid access information, and each grid has a unique code. All grids have similar shapes, and the lower grids are obtained by splitting the upper grids. (2) Divide the ephemeris into multiple segments, and determine whether the period meets the lighting conditions required for optical remote sensing satellite observation. The period that does not meet the lighting conditions is directly removed and no longer participates in the calculation; the specific method is: (201) The time period of the segmented ephemeris is set to 2 minutes, and each segment of the ephemeris is a set of subsatellite points SEG = {p1, p2, ..., p i-1 ,p i ,p i+1 ,……,p 120 }, for any subsatellite point p in SEG i , if any of the following conditions exist: ； ； Then p i The points after that are separated from SEG, generating a new segmented ephemeris SEG1={p i ,p i+1 ,……,p 120 }; (202) Using the satellite pointing calculation method, calculate each segment of the ephemeris {p1, p2, ..., p n }, the specific method is: The earth surface equation in the J2000 coordinate system is: (Formula 1) Among them, a and b are the equatorial radius of the earth, a=b=6378137m; c is the polar radius of the Earth, c=a×(1-f); f is the earth's flattening, f=1 / 298.257223563; Satellite boresight vector in J2000 coordinate system The parametric equation of the straight line vector is: (Formula 2) Among them, x0, y0, and z0 are the spatial positions of the satellite in the J2000 coordinate system at a certain time, which are obtained from the ephemeris file; when the satellite side viewing angle is α, m, n, and p are the satellite line of sight vectors In the three components of the coordinate axis, the specific calculation process is as follows: Satellite boresight in orbital coordinate system ; (Formula 3) Satellite boresight in J2000 coordinate system =(m,n,p)= ; (Formula 4) in, The transformation matrix from the orbital system to the J2000 coordinate system. The elements of this transformation matrix are: (Formula 5) Among them, i=1,2,3 corresponds to the three components of each row vector in the transformation matrix C; C(1,i), C(2,j) and C(3,k) are the row vectors of the first, second and third rows of the matrix C respectively; Represents the position vector of the satellite in the J2000 coordinate system, obtained from the ephemeris file; Represents the velocity vector in the satellite J2000 coordinate system, obtained from the ephemeris file; Combine (Formula 1)-(Formula 4) and solve the equations to calculate the intersection of the satellite's visual axis and the earth's surface, that is, the ground pointing point of the satellite at a certain time and with a certain field of view angle; In each ephemeris segment, the satellite's coverage of the ground is a quadrilateral, and the vertices of the quadrilateral are (P1, P2, P3, P4), where: Substitute the satellite J2000 spatial position and the minimum side viewing angle range at the beginning of each ephemeris into (Formula 1)-(Formula 4) to calculate the P1 point; Substitute the satellite J2000 spatial position and the maximum side viewing angle range at the beginning of each ephemeris into (Formula 1)-(Formula 4) to calculate the P2 point; Substitute the satellite J2000 spatial position and minimum side viewing angle range at the last moment of each ephemeris into (Formula 1)-(Formula 4) to calculate the P3 point; Substitute the satellite J2000 spatial position and the maximum side viewing angle range at the last moment of each ephemeris into (Formula 1)-(Formula 4) to calculate the P4 point; (203) If the satellite observation has requirements for lighting conditions, then determine the lighting conditions of the segmented ephemeris coverage range, and calculate the solar altitude angles sun1, sun2, sun3, and sun4 of points P1, P2, P3, and P4 respectively. If the four solar altitude angles are all less than 0, the area does not meet the observation conditions and is directly discarded; if the four solar altitude angles are all greater than or equal to 0, the area fully meets the imaging conditions; if there are K solar altitude angles greater than 0, K = 1, 2, 3, then the area partially meets the imaging conditions; The calculation formula for the solar altitude angle is: (Formula 6) Where: θ — the solar altitude angle at the observation point, in degrees; — latitude of the observation site, in degrees; δ — solar declination, in degrees; τ — observation geotime angle, in degrees; The calculation formula for the solar declination δ is: δ = 23.45 × sin [(360 / 365) × (284 + n)] (Formula 7) Where n is the number of days in a year, calculated from the time in the ephemeris; The calculation formula of the observation time angle τ is: (Formula 8) in: st is the true solar time, that is, the local time required for observation. For every 15 degrees of longitude on the earth's surface, the local time differs by 1 hour; Δt represents the time difference between the local time of the observation requirement and the Coordinated Universal Time; lon is the longitude of the observation requirement; t is the corresponding ephemeris time; tz is the time zone number corresponding to the ephemeris time standard. If the ephemeris time standard is Coordinated Universal Time, tz=0; if the ephemeris time standard is Beijing Time, tz=8; (204) Save the access information of the time period into the segmented filtering information list: SegList={s i |i=1,2,…}, where s i Contains the start time t1, end time t2 of the time period, and whether the imaging conditions are met; (3) Set the time step and the side view step. For the time period that meets the lighting conditions, calculate the access information of each grid that the satellite can cover in each time period according to the two dimensions of time step and side view step. The access information includes the access time and the side view angle. If the coverage area of ​​a time period partially meets the imaging conditions, calculate whether each grid meets the imaging conditions, and save the access information of the grids that meet the imaging conditions into the database. The specific method is as follows: (301) Select the satellite field of view angle as the step length of the side viewing angle to ensure that the satellite pointing crosses one grid each time the side viewing angle steps. When the satellite is in a high latitude area, the distance in the longitude direction becomes smaller, the width of the grid becomes smaller, and the side viewing angle steps once. If the satellite pointing crosses multiple grids and these grids are all within the range of the satellite field of view angle Vα, then the side viewing angles of these grids are the same. (302) Set the time step length to the ratio of the satellite observation width Vw to the distance between adjacent sub-satellite points, that is, the time step length ts = Vw / 7.8; if the satellite observation width Vw is less than 7.8, set ts = 1. In this case, if the satellite crosses multiple grids during one time step, the access time of these grids is still considered to be the same; (303) Traverse the pre-screening time period set SegList that meets the conditions, for any time period s i 、 Satellite side viewing range (α min , α max ), and process as follows: ①Let t = t1, α = α min ; ② When t>t2, execute step (304); ③When α>α max , execute step ⑧; ④ Calculate the ground pointing point p of the satellite at time t and side angle α; ⑤If s i If the imaging conditions of the time period are fully met, then execute step ⑥; if they are partially met, then calculate the solar altitude angle of point p at time t; if the solar altitude angle does not meet the imaging conditions, then execute step ⑦; ⑥ Calculate the grid number Gid where point p is located according to the method of step (104), record the grid access information GAI, set GAI.gid=Gid, GAI.t=t, GAI.α=α, and save the grid access information to the grid access information list GAIList; ⑦Let α=α+Vα, and execute step ③; ⑧Let t=t+ts, and execute step ②; (304) For high-latitude areas and grid widths less than 7.8 km, the grid access information of GAIList is checked. When the grid column number in a row is crossed, the access information of the crossed grid is supplemented, and the access information of the crossed grid is equal to the access information of the grid in the previous column in the row; when the grid row number in a column is crossed, the access information of the crossed grid is supplemented, and the access information of the crossed grid is equal to the access information of the grid in the previous row in the column; (305) Saving the access information to the database; in order to improve data retrieval efficiency, deleting outdated data, that is, deleting data whose access time is earlier than the current time; (4) Calculate the code of the grid where the observation requirement is located based on the longitude and latitude of the observation requirement, and retrieve the satellite's access information to the grid from the database as the satellite's access information to the observation requirement.

2. The satellite access information service method based on global gridding according to claim 1, characterized in that: The specific method of step (1) is: (101) The grid level n is determined according to the ground width observed by satellite, where n ≥ 2. The world is divided into a multi-level grid system according to longitude and latitude. The first level division method is: Horizontally, from the equator to the north and south poles, each latitude difference of 6° is a row, 15 rows in each hemisphere, a total of 30 rows; vertically, from longitude 0° to the east and west, each longitude difference of 6° is a column, 30 columns in each hemisphere, a total of 60 columns; Divide the first-level grid into 4 equal parts to obtain the second-level grid, and divide the second-level grid into 4 equal parts to obtain the third-level grid, and so on, until the n-th level grid; (102) Construct a grid coding rule, including layer number, northern and southern hemisphere code, row number, eastern and western hemisphere number, and column number, where: The layer number is a positive integer, the first layer is 1, the second layer is 2, and so on; The codes for the northern and southern hemispheres are: N for the northern hemisphere, S for the southern hemisphere; The row number is a positive integer, starting from 1, and increasing by 1 from the equator to the second level. The maximum value of the row number is different for different levels. The maximum value of the first level is 15, and it doubles with each level. The Eastern and Western Hemispheres are numbered as follows: E for the Eastern Hemisphere and W for the Western Hemisphere; The column number is a positive integer, starting from 1, and increasing by 1 from the 0 degree meridian to both sides. The maximum value of the column number is different at different levels. The maximum value of the first level is 30, and it doubles at each level. (103) Calculate the longitude and latitude range of the grid according to the grid number; if the layer number of the grid is lid, the row number is row, and the column number is col, then the longitude and latitude range of the grid is calculated as follows: The latitude and longitude difference of this layer grid is Δ = 6 / 2 lid-1 ; Grid minimum longitude minLon = (col - 1) × Δ; Grid maximum longitude maxLon = col×Δ; Grid minimum latitude minLat = (row - 1) × Δ; Grid maximum latitude maxLat = row×Δ; If the code for the eastern and western hemispheres is W, let minLon = -minLon, maxLon = -maxLon; If the code for the northern and southern hemispheres is S, let minLat = -minLat, maxLat = -maxLat; (104) According to the longitude and latitude coordinates (lon, Lat) and the layer number lid, calculate the grid number where it is located: The longitude and latitude difference of this layer grid is Δ = 6 / 2 lid-1 ; Grid column number col = TopInt(lon / Δ) Grid row number row = TopInt(Lat / Δ) TopInt is the upward rounding function; If lon ≥ 0, the code for the eastern and western hemispheres is E, otherwise, the code for the eastern and western hemispheres is W; If Lat ≥ 0, the code for the northern and southern hemispheres is N, otherwise, the code for the northern and southern hemispheres is S; (105) A grid layer with a grid width at the equator that matches the satellite observation width is used as the calculation carrier. The first layer of grids has 60 grids per row. The width of the first layer of grids at the equator is 666 km. The width of each subsequent layer of grids at the equator is half of the width of the previous layer of grids at the equator and is rounded up. (106) A two-layer storage architecture is used to store grid access information. The first layer contains a table named SatGridLevel, which stores the applicable grid levels of different satellites and is used to quickly determine the grid level corresponding to the satellite during grid access calculation; The second layer contains multiple tables, which store the satellite's access information to the grid. Each satellite corresponds to a table, and the table name is: satellite code + AI.

3. The satellite access information service method based on global gridding according to claim 1, characterized in that: The specific method of step (4) is: (401) After receiving the observation requirement, according to the resolution requirement in the observation requirement, obtain the applicable satellite, determine which grid access information table to query next, and put it into the table set to be queried TabList; (402) According to the observation requirement location, the grid number Gid corresponding to the observation requirement is calculated in the manner of step (104); (403) According to the grid number Gid, each table in TabList is queried to obtain the access satellite, access time and access angle required for observation.

Images