A real-time satellite coverage calculation method and system based on an equal-area spherical grid, an electronic device, and a storage medium

Abstract

The application provides a real-time satellite coverage calculation method and system based on an equal-area spherical grid, electronic equipment and a storage medium, and relates to the technical field of space mission planning. The method comprises the following steps: acquiring satellite orbit attitude parameters, extracting a visual axis direction and simplifying a visual axis trajectory into a plurality of discrete line segments along a time axis; initializing a nested structure pixel index of a target area based on a hierarchical rule of an equal-area spherical grid HEALPix; inputting the line segment data and the pixel index into a GPU, and decoding the pixel index into a direction vector in the GPU; calculating the intersection state of each pixel vector and a corresponding line segment in parallel in the GPU according to the satellite payload field of view shape, accumulating the coverage time; and outputting the coverage rate of the target area based on the coverage statistical result. The application eliminates high-latitude area distortion by using an equal-area grid, greatly reduces the calculation complexity by using a line segment trajectory approximation and GPU parallelism, and realizes high-precision calculation of a second-level real-time coverage rate of a user-drawn polygon.

Description

Technical Field

[0001] This invention relates to the field of spacecraft orbit control and mission planning technology, and in particular to a real-time satellite coverage calculation method, system, electronic device and storage medium based on an equal-area spherical grid. Background Technology

[0002] With the rapid development of satellite technology, satellite ground systems need to process massive amounts of Earth observation data. Among these, global or specific regional coverage statistics are core key indicators for evaluating the performance of communication and remote sensing satellite systems, mission planning, emergency observation, and online consultation. Users often need to define polygons of arbitrary shapes on Web-GIS systems and require the system to provide feedback within seconds on "how many times the area has been viewed by satellite and what the coverage rate is."

[0003] Currently, conventional statistical methods are mainly based on traditional latitude and longitude grids. For example, NASA's JPL open-source tool uses an equally spaced latitude and longitude grid (e.g., 0.1° × 0.1°) to convert the instantaneous field of view into polygons and uses a spatial database to determine whether a point is inside the polygon. However, this traditional technique has serious drawbacks: First, the area weighting is severely distorted. The area of ​​the latitude and longitude grid decreases sharply with latitude (for example, the area at 80°N is only 17% of that at the equator). In order to compensate for the area error, the system must introduce a cos(latitude) weight for correction. This coefficient has a nonlinear abrupt change at the edge of the polygon, resulting in the system error of high latitude coverage being 8% to 15% lower.

[0004] Secondly, the query and calculation efficiency is extremely low. The number of global grid points is enormous, and relying on database-based GIS calculations is not only slow in terms of querying, but also cannot handle concurrent users. Furthermore, if the ellipsoidal strip dynamic union method is used to calculate the coverage area, the time complexity of polygon union algorithms (such as Weiler-Atherton) is O(n^2). 2 When processing thousands of stripes a day, the computation time increases exponentially, often taking several minutes, and can even lead to memory overflow and explosion, making it impossible to meet the real-time requirements of online consultations such as "hand-drawn polygons and second-level response". Summary of the Invention

[0005] The purpose of this application is to overcome the dual defects of existing technologies, namely, "the distortion of high-latitude area statistics caused by latitude and longitude weight correction" and "the explosion of computational load caused by continuous union of polygons". This application proposes a real-time satellite coverage calculation method, system, electronic device and storage medium based on equal-area spherical grid. This application can realize the online closed-loop operation of "user arbitrarily hand-drawing polygons → obtaining extremely low error coverage in less than 1 second → driving task rearrangement".

[0006] The first aspect of this invention discloses a real-time satellite coverage calculation method based on an equal-area spherical grid, the method comprising: Step S1: Obtain satellite orbit and payload parameters, calculate satellite attitude and line-of-sight trajectory within a time period, and simplify the continuous line-of-sight trajectory into multiple discrete time segment data based on an error threshold along the time axis; Step S2: Obtain the target region. Based on the hierarchical partitioning rules of the equal-area spherical mesh HEALPix, initialize the nested structure pixel index corresponding to the target region, and transfer the nested structure pixel index and the time line segment data to the video memory of the graphics processor GPU. Step S3: Perform parallel computing within the GPU: decode the corresponding spherical direction vector based on the pixel index, establish a geometric model based on the shape of the satellite payload field of view, calculate the intersection or distance between each pixel direction vector and the instantaneous field of view corresponding to each time segment in parallel, and count the coverage status and cumulative coverage duration of each pixel. Step S4: Return the pixel coverage result output by the GPU to the central processing unit (CPU), calculate the overall coverage of the target area, and output it.

[0007] Furthermore, in step S1, the continuous line-of-sight trajectory is simplified along the time axis into multiple discrete time segment data, specifically including: The satellite position and attitude matrix in the Earth-centered Earth-fixed system is calculated using an orbit extrapolation model, and the line-of-sight direction vector of the sensor is extracted for each calculation step. A trajectory simplification algorithm based on angle thresholds is adopted to simplify the line segment of the line segment of the line axis direction vector sequence. When the maximum angular deviation between each sampling point in a trajectory and the line connecting the start and end points is less than a set threshold, the trajectory segment is approximated by the line segment represented by the start and end points, generating time line segment data containing start time, end time, start direction and end direction.

[0008] Furthermore, in step S2, the single-pixel area of ​​the equal-area spherical mesh HEALPix is ​​4πR² / (12×NSIDE²), where R is the Earth's radius and NSIDE is the mesh resolution parameter raised to the power of 2; the nested structure pixel index is in one-dimensional integer form, and there is no need to transmit three-dimensional coordinate data to the GPU.

[0009] Furthermore, when the satellite payload's field of view is a conical field of view, the time segment data includes the normalized direction vectors of the visual axes at both ends of the segment and the half-cone angle; the specific calculation process of step S3 includes: In the GPU worker thread, the nested pixel index is decoded into a direction vector V through bitwise operations; Traverse all time segments and perform a fast rejection test: calculate the dot product of the direction vector V and the view axes at both ends of the segment. If the dot product results are all less than the cosine of the half-cone angle, the segment is determined to be uncovered; if they are all greater than the cosine, the entire segment is determined to be covered; otherwise, treat the time segment as a great circle arc on a sphere, analytically solve the time sub-intervals where the distance between the direction vector V and the great circle arc is less than the half-cone angle, and accumulate them into the coverage duration of the pixel.

[0010] Furthermore, when the satellite payload's field of view is rectangular, the time segment data includes the direction vectors of the four corner points of the start and end points; the specific calculation process of step S3 includes: In the GPU worker thread, the spherical quadrilateral field of view is reconstructed for any interpolation moment corresponding to the time segment; A half-space test is performed on the four large circular edges of the spherical quadrilateral field of view. The cross product of adjacent corner points is calculated as the normal vector. If the direction vector V has the same sign as the dot product of the normal vectors of the four large circular edges, the pixel is determined to be inside the rectangular field of view. Sampling interpolation or binary search is performed in combination with the start and end times of the line segment to determine and accumulate the coverage time of the pixel.

[0011] Furthermore, before performing a half-space test on the rectangular field of view, a circumscribed cone containing the rectangular field of view is constructed in the GPU. A fast rejection test of the circumscribed cone is first performed on the direction vector V, and only the pixels that pass the test are subjected to a precise four-circle half-space test.

[0012] Furthermore, in step S4, when calculating the overall coverage of the target area, multi-resolution aggregation statistics and time-weighted statistics are supported; the instantaneous coverage is defined as the ratio of the number of covered pixels to the total number of pixels in the target area; the time-weighted coverage is defined as the sum of the products of the duration of each observation and the pixel coverage divided by the total observation time.

[0013] A second aspect of this invention discloses a real-time satellite coverage calculation system based on an equal-area spherical grid, the system comprising: The data preprocessing and segmentation module is used to acquire satellite orbit and payload parameters, calculate the line-of-sight trajectory within a time period, and simplify it into multiple discrete time segment data. The equal-area mesh initialization module is used to initialize the nested structure pixel index of the target region based on the equal-area spherical mesh HEALPix. The GPU parallel acceleration module is used to receive line segment data and pixel indexes, decode direction vectors in parallel, and calculate the coverage status and cumulative coverage time of each pixel in combination with the field shape. The coverage evaluation module receives GPU output results, calculates the overall coverage and related statistical indicators of the target area, and provides visualization output.

[0014] A third aspect of this invention discloses an electronic device. The electronic device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements the steps of the real-time satellite coverage calculation method based on an equal-area spherical grid according to any one of the first aspects of this disclosure.

[0015] A fourth aspect of this invention discloses a computer-readable storage medium. The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of a real-time satellite coverage calculation method based on an equal-area spherical grid, as described in any of the first aspects of this disclosure.

[0016] The beneficial effects of this application are as follows: 1. Revolutionary Improvement in Accuracy (Completely Eliminating Latitude Error): This application is the first to introduce HEALPix equal-area pixels from the astronomical field into satellite Earth coverage calculations, perfectly replacing latitude and longitude dot matrix with equal-area pixels. Regardless of whether the pixel is located at the equator or the poles, its area is strictly equal, thus completely eliminating the computationally intensive and easily distorted cos(θ) weight correction step. Measured data show that in the high-latitude region of the Arctic shipping route (76°~82°N), the coverage area error drops sharply from 11.4% using traditional methods to below 0.9%.

[0017] 2. Order-of-magnitude leap in response speed (transforming asynchronous into real-time): Traditional solutions feed the instantaneous field of view per second into the GPU for judgment, or perform polygon union on the CPU (taking nearly 5 minutes). This application cleverly utilizes the continuous and smooth physical characteristics of satellite orbital motion, approximating the entire trajectory with a small number of line segment polygonal envelopes (condensing 86,400 seconds of a day into thousands of line segments), and pushing the lightweight line segment data into the GPU. The computation time for a single high-concurrency query is reduced from hundreds of seconds to less than 0.1 seconds, fully supporting the "drag and drop to refresh" interactive experience of more than 50 users concurrently using the web interface. Attached Figure Description

[0018] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0019] Figure 1 A flowchart illustrating a real-time satellite coverage calculation method based on an equal-area spherical grid according to an embodiment of the present invention; Figure 2 This is an overall flowchart of a method provided in one embodiment of this application; Figure 3This is a structural diagram of a real-time satellite coverage calculation system based on an equal-area spherical grid according to an embodiment of the present invention; Figure 4 This is a structural diagram of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of the embodiments of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0021] The core of this invention lies in the fact that traditional calculations involving sampling points per second and finding intersections with polygons involve massive redundancy, while polygon Boolean unions are prone to causing dimensionality explosion. Since satellite orbital motion is relatively smooth, the movement of the field of view on the ground is a continuous process. Therefore, this application proposes a three-in-one solution: "HEALPix equal-area mesh + time-axis trajectory segmentation preprocessing + GPU fully parallel field of view intersection calculation".

[0022] The first aspect of this invention discloses a real-time satellite coverage calculation method based on an equal-area spherical grid. Figure 1 This is a flowchart illustrating a real-time satellite coverage calculation method based on an equal-area spherical grid according to an embodiment of the present invention. Figure 2 This is a flowchart illustrating the overall process architecture of a method provided in one embodiment of this application, as follows: Figure 1 and Figure 2 As shown, the method includes: Step S1: Obtain satellite orbit and payload parameters, calculate satellite attitude and line-of-sight trajectory within a time period, and simplify the continuous line-of-sight trajectory into multiple discrete time segment data based on an error threshold along the time axis; Specifically, the continuous line-of-sight trajectory is simplified along the time axis into multiple discrete time segment data, including: The satellite position and attitude matrix in the Earth-centered Earth-fixed system is calculated using an orbit extrapolation model, and the line-of-sight direction vector of the sensor is extracted for each calculation step. A trajectory simplification algorithm based on angle thresholds is adopted to simplify the line segment of the line segment of the line axis direction vector sequence. When the maximum angular deviation between each sampling point in a trajectory and the line connecting the start and end points is less than a set threshold, the trajectory segment is approximated by the line segment represented by the start and end points, generating time line segment data containing start time, end time, start direction and end direction.

[0023] A satellite day has 86,400 seconds; if we were to analyze it second by second, the computational workload would exceed hundreds of billions of times. The core of this step is to discretize the time: First, input the satellite orbital elements (TLE) or precise ephemeris. On the CPU, using the SGP4 analytical model or numerical integration algorithm, calculate the satellite position and velocity in the Earth-centered Earth-fixed frame (ECEF) with a step size of 1 second, and calculate the attitude quaternions by combining attitude modes such as Earth orientation.

[0024] Next, extract the optical axis (Z-axis) direction in the sensor's own system and transform it to the ECEF system. If it is a conical field of view (such as a communication antenna), record the direction of the optical axis and the half-cone angle; if it is a rectangular field of view (such as an optical remote sensing camera), extract the vectors of the four corner points.

[0025] Finally, the crucial line segment simplification is performed. The Ramer-Douglas-Peucker (RDP) algorithm or a greedy algorithm is used to simplify the continuous trajectory. If the center point deviation distance (maximum angular deviation) within a trajectory segment is less than a preset threshold (e.g., 0.1 degrees), then this trajectory spanning tens of seconds is replaced by a set of start-end "segments". After this step, the 86,400 seconds of trajectory is highly condensed into approximately 4,000 to 5,000 line segments, significantly reducing the size of the input data.

[0026] Step S2: Obtain the target region. Based on the hierarchical partitioning rules of the HEALPix equal-area spherical mesh, initialize the nested structure pixel index corresponding to the target region, and transfer the nested structure pixel index and the time segment data to the video memory of the graphics processing unit (GPU). The single pixel area of ​​the HEALPix equal-area spherical mesh is 4πR. 2 / (12×NSIDE 2 ), where R is the Earth's radius, NSIDE is the grid resolution parameter and is a power of 2; the nested structure pixel index is in one-dimensional integer form and there is no need to transmit three-dimensional coordinate data to the GPU.

[0027] Traditional 0.1° latitude and longitude grids exhibit significant area differences between high and low latitudes. This application employs the HEALPix (Hierarchical Equal Area isoLatitude Pixelation) algorithm. HEALPix is ​​based on dividing the celestial sphere into 12 equal-area base pixels. The total number of pixels N is determined by the resolution parameter NSIDE (e.g., 1024, i.e., the 10th layer). PIX =12×NSIDE 2 The key feature is that the area of ​​each pixel is strictly constant. .

[0028] For area coverage, instead of transmitting the 3D coordinate system of each grid, nested pixel indices (4-8 bytes of integers) are sent to the GPU memory. This significantly saves bus transmission bandwidth (accounting for only one-third of the floating-point coordinate data). In some specific embodiments, the specific implementation steps for initializing based on HEALPix equal-area meshes are as follows: Step S21, Mesh Generation Principle HEALPix divides the sphere into pixel units of equal area and equal intervals of latitude circles. Its resolution is controlled by the parameter NSIDE (values ​​are powers of 2, such as NSIDE = 2^k, k≥0). The formula for the total number of pixels is: NPIX = 12 × NSIDE².

[0029] The area of ​​each pixel is strictly equal, and its value is:

[0030] This feature eliminates the problem of area distortion in high-latitude regions caused by traditional latitude and longitude grids.

[0031] Resolution Selection: Determine the appropriate NSIDE value based on the spatial resolution and coverage area of ​​the satellite payload. For example, for satellite payloads with a large coverage area, NSIDE=256 (corresponding to approximately 780,000 pixels) and a pixel area of ​​approximately 53 km² can be selected; for high-resolution observations, NSIDE=1024 or higher can be selected.

[0032] The partitioning process uses a quadtree structure, and the specific steps are as follows: Basic division: The celestial sphere is divided into 12 base pixels of equal area using three mutually perpendicular large circles (corresponding to NSIDE=1).

[0033] Recursive subdivision: Each base pixel is continuously divided into four equal parts (doubling the NSIDE value each time) to form a hierarchical mesh. For example: When NSIDE=1, the number of pixels is 12; When NSIDE=2, the number of pixels is 12×4=48; When NSIDE=4, the number of pixels is 12×16=192; And so on.

[0034] During subdivision, a new great circle is formed by connecting the midpoint of the edge of each pixel to the center of the sphere, ensuring that the area of ​​all pixels remains equal after each subdivision. Sorting Method Selection: HEALPix offers two pixel sorting methods—RING and NESTED. RING sorting arranges pixels in latitudinal circles, suitable for spherical harmonic analysis; NESTED sorting has a hierarchical structure, suitable for multi-resolution analysis and fast neighborhood search. NESTED sorting is recommended for coverage statistics, facilitating detailed analysis of local areas.

[0035] Step S22: Coordinate Mapping and Pixel Index Calculation Let the latitude and longitude of the Earth's surface point be ( )( Longitude (latitude), needs to be converted to co-latitude in spherical coordinates. and longitude :

[0036] Calculate the pixel index using the ang2pix function of HEALPix. :

[0037] The pix2ang function can be used to reverse-map pixel indices to a spherical coordinate range:

[0038] Each pixel's boundary is enclosed by four great circle arcs, and its latitude and longitude range can be determined by interpolation.

[0039] Step S23: GPU initializes the mesh HEALPix is ​​a spherical discretization scheme with equal area. The nth layer has a partition of NSIDE = 2^n, and the total number of pixels N_pix = 12 × NSIDE². The 10th layer has NSIDE = 1024, a pixel angular resolution of approximately 3.4 angular minutes, and a ground resolution of approximately 3-10 kilometers (depending on orbital altitude), which is sufficient to meet the needs of communication coverage and remote sensing mission planning.

[0040] Pixels are stored using NESTED indices, rather than directly transmitting 3D coordinates. The indices are 32-bit or 64-bit integers (4-8 bytes), while the coordinates require 3 floating-point numbers (12 bytes). The GPU converts the indices into direction vectors in real time using an analytical formula, which is far less computationally expensive than the memory bandwidth savings, resulting in higher overall efficiency.

[0041] Step S3: Perform parallel computing within the GPU: decode the corresponding spherical direction vector based on the pixel index, establish a geometric model based on the shape of the satellite payload field of view, calculate the intersection or distance between each pixel direction vector and the instantaneous field of view corresponding to each time segment in parallel, and count the coverage status and cumulative coverage duration of each pixel. When the satellite payload's field of view is conical, the time segment data includes the normalized direction vectors and half-cone angles of the visual axes at both ends of the segment; the specific calculation process of step S3 includes: In the GPU worker thread, the nested pixel index is decoded into a direction vector V through bitwise operations; Traverse all time segments and perform a fast rejection test: calculate the dot product of the direction vector V and the view axes at both ends of the segment. If the dot product results are all less than the cosine of the half-cone angle, the segment is determined to be uncovered; if they are all greater than the cosine, the entire segment is determined to be covered; otherwise, treat the time segment as a great circle arc on a sphere, analytically solve the time sub-intervals where the distance between the direction vector V and the great circle arc is less than the half-cone angle, and accumulate them into the coverage duration of the pixel.

[0042] When the satellite payload's field of view is rectangular, the time segment data includes the direction vectors of the four corner points of the start and end points; the specific calculation process of step S3 includes: In the GPU worker thread, the spherical quadrilateral field of view is reconstructed for any interpolation moment corresponding to the time segment; A half-space test is performed on the four large circular edges of the spherical quadrilateral field of view. The cross product of adjacent corner points is calculated as the normal vector. If the direction vector V has the same sign as the dot product of the normal vectors of the four large circular edges, the pixel is determined to be inside the rectangular field of view. Sampling interpolation or binary search is performed in combination with the start and end times of the line segment to determine and accumulate the coverage time of the pixel.

[0043] Before performing a half-space test on the rectangular field of view, a circumscribed cone containing the rectangular field of view is constructed in the GPU. A fast rejection test of the circumscribed cone is first performed on the direction vector V, and a precise four-circle half-space test is performed only on the pixels that pass the test.

[0044] Specifically, after the GPU receives nested pixel indices and thousands of line segment data, it initializes a large number of concurrent threads, each processing one or a group of HEALPix pixels. Within the GPU registers, the integer indices are directly decoded into the corresponding direction vector V through bitwise operations and spherical mapping (eliminating the bottleneck of video memory access).

[0045] Based on different load fields of view, the GPU performs differentiated computations: Scenario A: Conical field of view (communication payload).

[0046] A line segment is defined as a great circle arc that moves at a constant speed from the starting point C0 to C1 along the view axis. For the current pixel vector V, the GPU executes: Fast Rejection Test: Calculate V C0 and V C1. If all values ​​are much smaller than cos(α) (α is the half-cone angle), skip the segment directly (no coverage); if all values ​​are much larger, the entire segment period is fully covered, and the time is added directly.

[0047] Analytical calculation: For line segments in the boundary state, construct a quadratic function approximation equation with respect to the interpolation time parameter λ, or directly use the spherical cosine theorem to find the minimum angular distance from the direction vector V to the great circle arc C0→C1. If the distance is less than the half-cone angle, calculate the critical time for entering and exiting the coverage and accumulate the time.

[0048] In some specific embodiments, for a conical field of view, the specific implementation is as follows: The geometric characteristics of a conical field of view are: the ground footprint is a spherical cap region with a circular boundary, uniquely determined by the view axis direction and the semi-cone angle. Coverage determination within a line segment is transformed into a spherical geometry problem: given a pixel direction vector V, the view axis along the line segment from C... Move to C Find when the angular distance between V and the motion line of sight is less than the half-cone angle α.

[0049] Step S311: Geometric Modeling and Analytical Solution The motion of the line of sight on the sphere is approximately a uniform great circle motion. For any time t∈[t... ,t The interpolation view axis direction is: C(t) = normalize( (1-λ)·C + λ·C ), Where λ = (tt) ) / (t -t ), The condition for pixel V to be covered is dot(V, C(t)). cos(α).

[0050] Since C(t) is a normalized linear interpolation, dot(V, C(t)) is a rational function of λ. To simplify the calculation, the following approximation can be used: Under the assumption of short line segments (the central angle corresponding to the line segment...) (5 degrees), the motion of the line of sight is approximated as a uniform linear motion in the plane, and the dot product is approximated as a quadratic function of λ: f(λ) = a·λ² + b·λ + c The coefficients are determined by the endpoint conditions. Solving f(λ) = cos(α) yields the critical λ values ​​for entering and exiting coverage, and thus the coverage time interval.

[0051] A more precise solution is to directly solve for the spherical triangle: consider the relationship between V and C. C Given a great circle arc, calculate the minimum angular distance from V to the arc. If the minimum distance is less than α, then a cover exists; the endpoints of the covered subinterval can be analytically determined using the sine and cosine laws for spheres.

[0052] Step S312, GPU parallel execution mode The compute shader takes the filtered pixel index as input, and each thread processes one pixel. The threads are organized into a one-dimensional or two-dimensional grid, and the worker thread group size is typically 256 or 512 to fully utilize the GPU's SIMT execution efficiency.

[0053] The thread execution flow is as follows: First, the NESTED index is decoded into spherical coordinates (θ, φ), and then converted into a three-dimensional unit vector V. Decoding utilizes the nested structure of HEALPix, extracting the face number (F=0-11) and pixel coordinates (x, y) through bitwise operations, and then obtaining the direction using the projection formula. This calculation is entirely performed within registers, requiring no access to video memory.

[0054] Then, initialize the cumulative coverage time of the pixel to zero. Traverse all line segments (approximately 5000, loaded from read-only memory or texture cache): For each line segment, perform a fast rejection test—calculate the dot product of V and the view axes at both ends of the line segment. If both are less than cos(α), the line segment has no coverage; if both are greater than cos(α), the entire segment is covered; otherwise, analytically solve for the coverage sub-intervals and accumulate the duration.

[0055] Finally, the cumulative coverage time (or coverage = time / 86400) is written to the output buffer. If only statistics (such as global average coverage) are needed, parallel reduction can be performed within the GPU to avoid large amounts of data backhaul.

[0056] Step S313, Accuracy and Error Control The error in approximating a line segment stems from replacing a smooth line-of-sight trajectory with a straight line segment. This error is related to the line segment length and orbital curvature: for low-Earth orbit satellites (altitude 500km), the orbital period is approximately 90 minutes, 20 seconds corresponds to an orbital arc length of approximately 140km, and the line-of-sight direction changes by about 1-2 degrees. The introduced coverage boundary error is less than the pixel size (3-10km), satisfying the requirements... 1% area accuracy requirement.

[0057] For higher precision, adaptive encryption can be used: automatically shorten the line segment length during orbital maneuvers or periods of rapid attitude change; allow for longer line segments during periods of stable flight. Alternatively, secondary subdivision verification can be performed on boundary pixels (covering time near the endpoints of the line segments).

[0058] Scene B: Rectangular field of view (remote sensing payload).

[0059] The geometry of a quadrilateral is more complex; each line segment contains three-dimensional vectors of its eight corner points, from the start to the end. GPU execution: Early rejection filtering: Construct a rough circumscribed cone containing quadrilaterals, and exclude anything not within this cone.

[0060] Spherical Quadrilateral Reconstruction and Half-Space Testing: For candidate pixels, a spherical polygon is reconstructed at the sampling times of the line segment (e.g., start point, midpoint, end point). The cross product of adjacent corner points yields the normal vectors of the four boundary great circles. If the dot product of the pixel vector V and these four normal vectors has the same sign, it means the pixel is completely within the field of view. If the state is reversed at any given time, a binary search is used to quickly approximate the critical moment when the pixel actually passes through the field of view.

[0061] In some specific embodiments, for a rectangular field of view, the specific implementation is as follows: The geometric complexity of a rectangular field of view lies in the fact that the ground footprint is a spherical quadrilateral, the boundary of which is formed by connecting the four corner points with great circle arcs, and its shape and orientation may change over time. This requires richer line segment representations and more complex coverage determinations.

[0062] Step S321, Line Segment Data Structure Each line segment needs to describe the complete motion of the rectangle from the start time to the end time. It uses a four-corner point independent trajectory storage method, storing the directions of the four corner points of the start and end points (a total of 8 three-dimensional vectors). Within the GPU, for any time t, the four corner points are interpolated to reconstruct the spherical quadrilateral. The data size is 8 × 3 × 4 = 96 bytes / line segment, approximately 480KB for 5000 segments.

[0063] Step S322: Rectangular Coverage Determination and Solution For a fixed pixel V and time t, the reconstructed rectangle is composed of four corner points P. (t),P (t),P (t),P (t) Definition. Coverage determination is transformed into testing the position of a spherical point within a spherical polygon.

[0064] The standard method is the half-space test: For each of the four sides of the rectangle, there is a great circle whose normal vector is the cross product (normalized) of the adjacent corner points. If the dot product of V and the normal vectors of the four sides has the same sign (both positive or all negative, depending on the order of the vertices), then V is inside the rectangle.

[0065] For continuous coverage within a line segment, we need to solve for the times when V enters and exits the rectangle. Since the rectangle boundaries change over time, analytical solutions are quite complex. Practical strategies include: Sampling-interpolation method: Reconstruct and test the rectangle at three time points: the start point, midpoint, and end point of the line segment. If the coverage state remains unchanged (all present or all absent), make a direct determination; if it changes, subdivide the area near the change (e.g., binary search) to sufficient precision.

[0066] Boundary trajectory method: Track the movement of the four sides of a rectangle on a sphere, transforming the problem into the intersection detection of V with the moving great circle. Each side is solved independently, and then the solutions are combined to determine the coverage area.

[0067] In GPU implementations, it is recommended to combine sampling-interpolation with early rejection: first use an circumscribed cone for quick screening, and then perform precise rectangle tests only on pixels that may be covered, reducing computational load.

[0068] Step S4: Return the pixel coverage result output by the GPU to the central processing unit (CPU), calculate the overall coverage of the target area, and output it.

[0069] In step S4, when calculating the overall coverage of the target area, multi-resolution aggregation statistics and time-weighted statistics are supported; the instantaneous coverage is defined as the ratio of the number of covered pixels to the total number of pixels in the target area; the time-weighted coverage is defined as the sum of the products of the duration of each observation and the pixel coverage, divided by the total observation time.

[0070] In some specific embodiments, the specific process for calculating coverage statistics is as follows: Step S411: Calculation of instantaneous coverage in a single instance Let the spherical area covered by a single satellite observation be... The set of HEALPix pixels it covers is:

[0071] Instantaneous coverage is defined as the ratio of the number of covered pixels to the total number of pixels:

[0072] Since the pixel areas are equal, this ratio is an unbiased estimate of the area coverage. Step S412: Calculation of cumulative coverage over time Within the time interval [t0, tn], suppose the satellite performs m observations, with each observation timestamp being... The set of covered pixels is ( Cumulative coverage calculation is divided into: Cumulative coverage pixel set: Cumulative coverage:

[0073] Step S413, Time-weighted coverage If we consider the duration Δtj of each observation, then the time-weighted coverage is:

[0074] The total observation time is: .

[0075] In some specific embodiments, the specific implementation methods of multi-resolution analysis and adaptive statistics are as follows: Step S421, Multi-resolution aggregation HEALPix's hierarchical structure supports fast aggregation statistics. For example, a single pixel at low resolution (NSIDE=8) can be divided into four sub-pixels at high resolution (NSIDE=16). Coverage can be aggregated from high resolution to low resolution.

[0076] in For low-resolution pixels, Its sub-pixels.

[0077] Step S422, Regional Adaptive Statistics For a specific region R (such as the Arctic region), the pixel set PR within it is statistically analyzed, and the region coverage is:

[0078] Example: If the Arctic region contains 120 pixels (NSIDE=32), and 80 of them are covered, then CR=80 / 120=66.7%.

[0079] In some specific embodiments, the specific implementation methods of error analysis and accuracy control are as follows: Step S431: Relationship between resolution selection and error The statistical error mainly originates from pixel discretization. Let the diameter of the satellite payload footprint be D (radians), then the approximate side length of a HEALPix pixel is: ; To reduce edge effects, the following must be met: General requirements: ; Example: If D=1 If the value is approximately 0.0175 radians, then NSIDE must be ≥90 (take the nearest power of 2, such as NSIDE=128).

[0080] Step S432, Boundary pixel processing When the observation area only partially covers pixels, weighted allocation can be used: Binary method: If the center of a pixel is within the observation area, it is counted as 1; otherwise, it is counted as 0 (simple but may overestimate the error). Area ratio method: Assign weights w∈[0,1] according to the overlapping area ratio, and calculate the coverage rate using weighted statistics. In some specific embodiments, the visualization and output are implemented as follows: Step S441, Projection Conversion HEALPix data is often displayed using the Mollweide projection, a constant-area projection that preserves the pixel area and avoids visual distortion. The projection formula is:

[0081] in satisfy .

[0082] Step S442: Output of statistical results The final global coverage matrix C(ipix) is generated, where each pixel stores the number of times it is covered or the percentage of time it is covered, and can be exported in FITS or NetCDF format.

[0083] To verify the advancement of this invention, experiments were conducted on a server equipped with an 8-core CPU and a consumer-grade RTX GPU. The experiment targeted a 2000 square kilometer Arctic observation target at 76°–82°N: (1) Using the well-known technique A (Latitude and Longitude PostgreSQL): It takes 2.3 seconds. Due to the extremely dense edge points and severe cosine decay, the statistical coverage has an absolute error of 11.4%.

[0084] (2) Using the well-known technique B (elliptic strip union method): the error is 0.3%, but the execution time for a single concurrent operation is as high as 285 seconds.

[0085] (3) The proposed solution has an error of only 0.9% and an execution time of only 0.08 seconds. It completely fills the technological gap between "high-precision polar statistics" and "interactive second-level response".

[0086] In summary, this application achieves accurate calculation of global coverage through HEALPix's equal-area grid generation, coordinate mapping, multi-resolution statistics, and error control. This method ensures fairness in high- and low-latitude statistics and is suitable for satellite payload performance evaluation and mission planning.

[0087] The second aspect of this invention discloses a real-time satellite coverage calculation system based on an equal-area spherical grid. Figure 3 This is a structural diagram of a real-time satellite coverage calculation system based on an equal-area spherical grid according to an embodiment of the present invention; as shown. Figure 3 As shown, the system 100 includes: The data preprocessing and segmentation module 101 is used to acquire satellite orbit and payload parameters, calculate the line-of-sight trajectory within a time period, and simplify it into multiple discrete time segment data. Specifically, it is responsible for receiving the ephemeris, time range, and latitude-longitude polygon input by the user, extrapolating the trajectory using the SGP4 model and attitude matrix, and packaging it into thousands of line segment structures using a threshold smoothing algorithm. The equal-area mesh initialization module 102 is used to initialize the nested structure pixel index of the target region based on the equal-area spherical mesh HEALPix; specifically, it determines the NSIDE level of HEALPix according to the system accuracy requirements, delineates all NESTED pixel numbers involved in the polygon, and allocates them as a one-dimensional array. The GPU parallel acceleration module 103 is used to receive line segment data and pixel index, decode direction vectors in parallel, and calculate the coverage status and cumulative coverage time of each pixel in combination with the field shape. It can concurrently decode vectors in the Compute Shader layer through memory scheduling (such as WebGPU writing), and quickly count intersections using the early rejection algorithm based on the algebraic rules of dot product and cross product. The coverage comprehensive evaluation module 104 is used to receive the GPU output results, calculate the overall coverage and related statistical indicators of the target area and output them in a visual format; perform multi-resolution aggregation on the dense one-dimensional array returned by the GPU (such as folding high-resolution pixel data into low-resolution areas for macroscopic display), and output a high-precision report free from weight interference.

[0088] A third aspect of this invention discloses an electronic device. The electronic device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements the steps of the real-time satellite coverage calculation method based on an equal-area spherical grid, as disclosed in any of the first aspects of this invention.

[0089] Figure 4 This is a structural diagram of an electronic device according to an embodiment of the present invention, such as... Figure 4As shown, the electronic device includes a processor, memory, communication interface, display screen, and input device connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, carrier networks, Near Field Communication (NFC), or other technologies. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, buttons, a trackball, or a touchpad mounted on the device's casing, or an external keyboard, touchpad, or mouse.

[0090] Those skilled in the art will understand that Figure 4 The structure shown is merely a structural diagram of the part related to the technical solution of this disclosure and does not constitute a limitation on the electronic device to which the solution of this application is applied. The specific electronic device may include more or fewer components than shown in the figure, or combine certain components, or have different component arrangements.

[0091] A fourth aspect of this invention discloses a computer-readable storage medium. The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of a real-time satellite coverage calculation method based on an equal-area spherical grid, as described in any of the first aspects of this invention.

[0092] Please note that the technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments have been described. However, as long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification. The above embodiments only illustrate several implementation methods of this application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims.

[0093] The above are preferred embodiments of the present invention. It should be noted that, for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A real-time satellite coverage calculation method based on an equal-area spherical grid, characterized in that, The method includes: Step S1: Obtain satellite orbit and payload parameters, calculate satellite attitude and line-of-sight trajectory within a time period, and simplify the continuous line-of-sight trajectory into multiple discrete time segment data based on an error threshold along the time axis; Step S2: Obtain the target region. Based on the hierarchical partitioning rules of the equal-area spherical mesh HEALPix, initialize the nested structure pixel index corresponding to the target region, and transfer the nested structure pixel index and the time line segment data to the video memory of the graphics processor GPU. Step S3: Perform parallel computing within the GPU: decode the corresponding spherical direction vector based on the pixel index, establish a geometric model based on the shape of the satellite payload field of view, calculate the intersection or distance between each pixel direction vector and the instantaneous field of view corresponding to each time segment in parallel, and count the coverage status and cumulative coverage duration of each pixel. Step S4: Return the pixel coverage result output by the GPU to the central processing unit (CPU), calculate the overall coverage of the target area, and output it.

2. The real-time satellite coverage calculation method based on an equal-area spherical grid according to claim 1, characterized in that, In step S1, the continuous line-of-sight trajectory is simplified along the time axis into multiple discrete time segment data, specifically including: The satellite position and attitude matrix in the Earth-centered Earth-fixed system is calculated using an orbit extrapolation model, and the line-of-sight direction vector of the sensor is extracted for each calculation step. A trajectory simplification algorithm based on angle thresholds is adopted to simplify the line segment of the line segment of the line axis direction vector sequence. When the maximum angular deviation between each sampling point in a trajectory and the line connecting the start and end points is less than a set threshold, the trajectory segment is approximated by the line segment represented by the start and end points, generating time line segment data containing start time, end time, start direction and end direction.

3. The real-time satellite coverage calculation method based on an equal-area spherical grid according to claim 2, characterized in that, In step S2, the area of ​​a single pixel of the equal-area spherical mesh HEALPix is: ， Where R is the Earth's radius, NSIDE is the grid resolution parameter and is a power of 2; the nested structure pixel index is in one-dimensional integer form and there is no need to transmit three-dimensional coordinate data to the GPU.

4. The real-time satellite coverage calculation method based on an equal-area spherical grid according to claim 1, characterized in that, When the satellite payload's field of view is conical, the time segment data includes the normalized direction vectors and half-cone angles of the visual axes at both ends of the segment; the specific calculation process of step S3 includes: In the GPU worker thread, the nested pixel index is decoded into a direction vector V through bitwise operations; Traverse all time segments and perform a fast rejection test: calculate the dot product of the direction vector V and the view axes at both ends of the segment. If the dot product results are all less than the cosine of the half-cone angle, the segment is determined to be uncovered; if they are all greater than the cosine of the half-cone angle, the entire segment is determined to be covered; otherwise, treat the time segment as a great circle arc on a sphere, analytically solve the time sub-intervals where the distance between the direction vector V and the great circle arc is less than the half-cone angle, and accumulate them into the coverage duration of the pixel.

5. The real-time satellite coverage calculation method based on an equal-area spherical grid according to claim 1, characterized in that, When the satellite payload's field of view is rectangular, the time segment data includes the direction vectors of the four corner points of the start and end points; the specific calculation process of step S3 includes: In the GPU worker thread, the spherical quadrilateral field of view is reconstructed for any interpolation moment corresponding to the time segment; A half-space test is performed on the four large circular edges of the spherical quadrilateral field of view. The cross product of adjacent corner points is calculated as the normal vector. If the direction vector V has the same sign as the dot product of the normal vectors of the four large circular edges, the pixel is determined to be inside the rectangular field of view. Sampling interpolation or binary search is performed in combination with the start and end times of the line segment to determine and accumulate the coverage time of the pixel.

6. The real-time satellite coverage calculation method based on an equal-area spherical grid according to claim 5, characterized in that, Before performing a half-space test on the rectangular field of view, a circumscribed cone containing the rectangular field of view is constructed in the GPU. A fast rejection test of the circumscribed cone is first performed on the direction vector V, and a precise four-circle half-space test is performed only on the pixels that pass the test.

7. The real-time satellite coverage calculation method based on an equal-area spherical grid according to claim 1, characterized in that, In step S4, when calculating the overall coverage of the target area, multi-resolution aggregation statistics and time-weighted statistics are supported; the instantaneous coverage is defined as the ratio of the number of covered pixels to the total number of pixels in the target area; the time-weighted coverage is defined as the sum of the products of the duration of each observation and the pixel coverage, divided by the total observation time.

8. A real-time satellite coverage calculation system based on an equal-area spherical grid, characterized in that, The system includes: The data preprocessing and segmentation module is used to acquire satellite orbit and payload parameters, calculate the line-of-sight trajectory within a time period, and simplify it into multiple discrete time segment data. The equal-area mesh initialization module is used to initialize the nested structure pixel index of the target region based on the equal-area spherical mesh HEALPix. The GPU parallel acceleration module is used to receive line segment data and pixel indexes, decode direction vectors in parallel, and calculate the coverage status and cumulative coverage time of each pixel in combination with the field shape. The coverage evaluation module receives GPU output results, calculates the overall coverage and related statistical indicators of the target area, and provides visualization output.

9. An electronic device, characterized in that, The electronic device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements the steps of the real-time satellite coverage calculation method based on an equal-area spherical grid as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of the real-time satellite coverage calculation method based on an equal-area spherical grid as described in any one of claims 1 to 7.

Images