Osgearth-based satellite orbit eci anchor rendering method
Patent Information
- Application Number
- CN202611076339.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2026-07-20
- Publication Date
- 2026-09-11
- Estimated Expiration
- 2046-07-20
AI Technical Summary
[0007]有鉴于此,本发明旨在提出基于osgEarth的卫星轨道ECI锚定渲染方法,以解决如下技术问题:(1)如何在osgEarth场景图中实现卫星轨道的无漂移渲染,避免每帧重新进行SGP4采样所带来的巨大CPU计算开销;(2)如何在GPU端高效消除高轨道卫星轨道线在地球背面产生的错误显示,解决远距离尺度下深度缓冲精度不足导致的穿透伪影问题;(3)如何将卫星轨道与地面实体统一到同一时空基准,解决天地实体时空基准不一致导致的对齐问题
(1)本发明通过建立ECI惯性系静态存储轨道几何与MatrixTransform节点每帧GMST旋转相结合的存一次、转每帧的管线架构,将轨道线漂移问题转化为固定矩阵更新操作。主线程每帧仅需执行一次旋转矩阵更新,计算开销与场景中的卫星数量无关,显著降低了CPU负载,实现了大规模卫星场景下的无漂移实时渲染。同时,通过动态确定采样点数,确保了不同周期轨道的几何渲染精度。
Smart Images

Figure CN122597740B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of computer graphics and 3D geographic information system visualization technology, and in particular relates to a satellite orbit ECI anchoring rendering method based on osgEarth. Background Technology
[0002] In the fields of space situation simulation and digital twins, it is necessary to render satellite orbit lines in real time within a 3D Earth scene and display them collaboratively with ground vehicles, aircraft, ships, and other entities within the same scene. With the development of computer graphics and 3D geographic information system visualization technologies, higher demands are placed on the joint simulation and visualization of large-scale space and ground targets. However, existing technologies have the following shortcomings in achieving these functions.
[0003] First, there is the issue of orbital drift. Satellite orbits are represented as fixed ellipses relative to stars in the ECI geocentric inertial coordinate system, while scene graphs in mainstream 3D Earth engines like osgEarth are typically rooted in the ECEF geocentric geofixed coordinate system. Existing solutions usually store orbital sampling points directly in the ECEF coordinate system. Due to the Earth's rotation, orbital lines sampled at fixed epochs will continuously drift with the Earth's rotation. To address the drift problem, some solutions choose to resample SGP4 orbital propagation data every frame to update the orbital geometry. However, for large-scale scenes containing dozens or even hundreds of satellites, performing tens of thousands of SGP4 propagations per frame would generate enormous CPU computational overhead, making real-time rendering impractical. Current technologies lack an efficient pipeline architecture that statically stores orbital geometry in the ECI coordinate system and achieves drift-free rendering by rotating scene graph transformation nodes every frame.
[0004] Secondly, there is the issue of handling Earth occlusion from high-orbit satellites. Geostationary orbit satellites orbit at an altitude of approximately 35,786 kilometers, and their orbital lines cross the Earth's surface in space. In 3D rendering, the standard depth buffer lacks sufficient accuracy at such large scales, making it difficult to correctly distinguish the spatial relationship between the orbital lines and the Earth's surface, thus producing visual artifacts of the orbital lines penetrating the Earth's surface. Existing solutions typically use CPU-side logic to handle Earth occlusion, failing to offload occlusion detection to the GPU graphics pipeline for parallel execution. Current technology lacks an efficient solution to directly perform Earth occlusion culling within the GPU fragment shader to address the insufficient accuracy of the depth buffer.
[0005] Finally, there is an inconsistency between the spacetime and ground time references. Satellite SGP4 orbital propagation typically uses Coordinated Universal Time (UTC) as input, and the output ECI coordinates require Greenwich Mean Sidereal Time (GMT) rotation to obtain ECEF coordinates. Ground entities, on the other hand, usually use geographic coordinates, which are converted to world coordinates via the osgEarth geographic transformation node. The inconsistency in the coordinate systems and time references upon which these two transformation links rely results in a lack of a unified simulation time axis driving mechanism in the system. Current technology lacks a mechanism to unify satellite orbital propagation and ground entity kinematic models to the same simulation time reference, making it difficult to guarantee accurate alignment of spacetime and ground entities under complex spacetime transformations.
[0006] Therefore, it is urgent to develop a satellite orbit ECI anchoring rendering method based on osgEarth to solve the aforementioned technical problems such as orbit drift, occlusion artifacts, and spatiotemporal inconsistency. Summary of the Invention
[0007] In view of this, the present invention aims to propose a satellite orbit ECI anchoring rendering method based on osgEarth to solve the following technical problems: (1) how to achieve drift-free rendering of satellite orbits in the osgEarth scene map, avoiding the huge CPU computing overhead caused by re-sampling SGP4 every frame; (2) how to efficiently eliminate the erroneous display of high-orbit satellite orbit lines on the back of the Earth on the GPU side, and solve the problem of penetration artifacts caused by insufficient depth buffer accuracy at long distance scales; (3) how to unify satellite orbits and ground entities to the same spatiotemporal reference, and solve the alignment problem caused by the inconsistency of spatiotemporal references between the ground and earth entities.
[0008] To achieve the above objectives, the technical solution of the present invention is implemented as follows: This invention provides a satellite orbit ECI anchoring rendering method based on osgEarth, comprising the following steps: Step 1, ECI Inertial Frame Static Storage and GMST Rotation Rendering: In the ECI coordinate system, sample the orbital geometry data for half an orbital period before and after the current simulation time, and store it as a point set P_eci. The number of sampling points N is dynamically determined according to the orbital period T_orb. Construct a scene graph node hierarchy structure, which consists of a MapNode (map root node), an independent parent node P_orbit, a MatrixTransform node M_orbit, and orbital line geometry nodes. For each frame, construct a rotation matrix based on the GMST angle calculated at the current simulation time, and map the entire ECI orbital point set to the ECEF world coordinate system. Only update the transformation matrix of the MatrixTransform node M_orbit for each frame. Step 2, Anchor point translation correction: Monitor the deviation between the actual ECI position of the satellite and the orbital geometric center anchor point p_anchor. When the deviation exceeds the recentering threshold, calculate the offset vector and perform an overall vector addition translation on all N points in the point set P_eci. Step 3, Dual-modal orbit display and asynchronous calculation: Supports switching between inertial orbit mode and nadir trajectory mode. Orbit geometry sampling is executed asynchronously in a background thread. Timing coordination between the main thread and the background thread is completed through asynchronous request identifier generation, transmission and matching filtering. Step 4, GPU fragment shader Earth occlusion culling: Turn off depth testing in the rendering pipeline, implement the intersection determination between the viewpoint and the Earth sphere in the custom GLSL fragment shader, calculate the intersection point of the ray from the viewpoint to the current fragment and the Earth sphere, and perform the discard operation when the fragment is behind the Earth. Step 5: Unify the simulation timeline and align it with the spacetime of Earth and Earth: Using the simulation reference time plus the simulation elapsed time as the unified time reference, the satellite SGP4 propagation and the ground entity kinematic model use the same simulation time drive, and all entities are unified to the Earth scene map through the geographic transformation node of osgEarth.
[0009] Furthermore, in step one, the number of sampling points N is dynamically determined based on the orbital period T_orb using the following formula: ; Wherein, K_samp is the sampling density coefficient, N_max is the upper limit of the number of sampling points, and N_min is the lower limit of the number of sampling points; sampling is performed symmetrically around the current simulation time, and the current ECI position of the satellite is used as the midpoint anchor point of the polyline and forcibly inserted into the middle position of the point set.
[0010] Furthermore, in step one, a rotation matrix is constructed for each frame based on the GMST angle calculated at the current simulation time, mapping the entire ECI orbit point set to the ECEF world coordinate system. Specifically, this includes: Construct the rotation matrix R_z(θ_gmst) based on the GMST angle θ_gmst calculated at the current simulation time: ; The transformation matrix M for mapping the entire ECI orbit point set to the ECEF world coordinate system is calculated as follows: ; Where T(p_ecef) is the translation matrix of the satellite's current ECEF position.
[0011] Furthermore, in step two, the deviation between the satellite's actual ECI position and the orbital geometric center anchor point p_anchor is monitored. When the deviation exceeds the recentering threshold, the offset vector is calculated, and a global vector addition translation is performed on all N points in the point set P_eci. Specifically, this includes: Calculate the deviation: ; When δ exceeds the recentering threshold ε, calculate the offset vector: ; Perform global vector addition and translation on all N points in the point set P_eci: ; Where p_i is the i-th sampling point before translation, and p_i' is the i-th sampling point after translation.
[0012] Furthermore, in step three, the timing coordination between the main thread and the background thread is completed through asynchronous request identifier generation, transmission, and matching filtering, specifically including: Each time a mode switch or forced refresh occurs, the main thread increments the asynchronous request identifier requestId and passes the requestId along with the current track display mode type to the background thread; After the background thread completes the sampling, it returns the requestId, track display mode type, and calculation result together. When consuming results, the main thread performs two checks in sequence: the first check checks whether the returned requestId is equal to the latest request identifier, and the second check checks whether the returned track display mode type is consistent with the current display mode; after both checks pass, the calculation results are applied to the track geometry nodes of the scene graph.
[0013] Furthermore, in step four, the intersection determination between the viewpoint and the Earth sphere is implemented in the custom GLSL fragment shader, specifically including: The shader takes the Earth's center c_earth and Earth's radius R_earth as input parameters and calculates the ray direction: ; And the projected length of the line of sight to the center of the Earth: ; When t_c>0, the shortest distance d_min from the line of sight to the center of the Earth satisfies: ; If d_min 2 <R_earth 2 Then calculate the intersection parameters of the points entering the sphere: ; When t_enter > 0 and t_enter < |po| - ε_t, a discard operation is performed; where ε_t is the tangent direction tolerance.
[0014] Furthermore, in step five, using the simulation reference time plus the simulation elapsed time as a unified time reference, the satellite SGP4 propagation and ground entity kinematic models use the same simulation time drive, specifically: Satellite passes by: ; SGP4 orbital propagation is performed, and the output ECI coordinates are rotated by GMST to obtain ECEF coordinates; the ground entity uses Δt_sim as the time step to drive the kinematic model and physics engine. In addition, when the camera altitude is below the atmospheric hiding threshold H_hide, the orbital lines are forcibly hidden.
[0015] Furthermore, in step three, the dual-mode orbital display includes: In inertial orbit mode, the orbit line appears as an ellipse fixed relative to inertial space, rotating as a whole with the Earth's rotation. In the nadir trajectory mode, the sampled geographic coordinates are converted from osgEarth to world coordinates and appear in a figure-eight or sine shape as the Earth rotates.
[0016] Compared with existing technologies, the satellite orbit ECI anchoring rendering method based on osgEarth of the present invention has the following advantages: (1) This invention transforms the orbital drift problem into a fixed matrix update operation by establishing a pipeline architecture that combines static storage of orbital geometry in the ECI inertial frame with GMST rotation of the MatrixTransform node per frame, storing once and rotating every frame. The main thread only needs to perform the rotation matrix update once per frame, and the computational overhead is independent of the number of satellites in the scene, significantly reducing the CPU load and realizing drift-free real-time rendering in large-scale satellite scenes. At the same time, by dynamically determining the number of sampling points, the geometric rendering accuracy of orbits with different periods is ensured.
[0017] (2) The present invention designs an anchor point translation correction mechanism, which corrects position deviations through batch vector addition within a small time scale, avoiding the computational overhead caused by frequently re-executing SGP4 propagation. This mechanism adopts a dual-path collaborative strategy of anchor point translation as the fast path and asynchronous full-circle resampling as the slow path, which effectively corrects the floating-point errors accumulated over long-term operation while ensuring real-time rendering performance, thus ensuring the long-term accuracy of orbit geometry.
[0018] (3) This invention proposes an Earth occlusion culling scheme based on GPU fragment shaders. Under the condition of disabling depth testing, the occlusion detection is offloaded to the GPU parallel pipeline for execution by geometric determination of the intersection of the line of sight and the Earth sphere. This scheme effectively overcomes the artifact problem of orbit lines penetrating the Earth's surface caused by insufficient accuracy of the standard depth buffer in high orbit scenes, and the occlusion judgment overhead is constant, only increasing the computational load of fragment shaders by a small amount, resulting in high rendering efficiency.
[0019] (4) This invention solves the problem of inconsistent spatiotemporal references between the satellite SGP4 orbit propagation and the ground entity kinematic model by constructing a unified simulation time axis and integrating them with the same time reference. By combining the matrix transformation and geographic transformation mechanism of the hierarchical combination in the scene graph nodes, it ensures that the spatial relationship between the satellite body and the orbit line remains consistent under complex spatiotemporal transformations, and realizes high-precision simulation and display of space-ground collaboration. Attached Figure Description
[0020] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 This is a flowchart of the satellite orbit ECI anchoring rendering method based on osgEarth provided in the embodiments of the present invention; Figure 2 This is a schematic diagram of the scene graph node hierarchy structure construction in the ECI anchoring mechanism provided in this embodiment of the invention; Figure 3 This is a schematic diagram of the dual-path coordination mechanism in the ECI anchoring mechanism provided in this embodiment of the invention; Figure 4 This is a comparative schematic diagram of the dual-modal track display provided in an embodiment of the present invention; Figure 5 This is a schematic diagram of asynchronous computing timing provided in an embodiment of the present invention; Figure 6 This is a schematic diagram illustrating the principle of Earth Occlusion Canceling in the GPU Fragment Shader provided in this embodiment of the invention; Figure 7 This is a schematic diagram of the simulation timeline in the space-ground linkage timing provided in the embodiment of the present invention; Figure 8 This is a schematic diagram of the satellite orbit layer in the space-ground linkage timing provided in the embodiments of the present invention; Figure 9 This is a schematic diagram of the ground entity layer in the time-series of the space-ground linkage provided in the embodiment of the present invention. Detailed Implementation
[0021] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.
[0022] In the description of this invention, it should be understood that the terms center, longitudinal, transverse, upper, lower, front, rear, left, right, vertical, horizontal, top, bottom, inner, and outer, indicating orientation or positional relationships, are based on the orientation or positional relationships shown in the accompanying drawings and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, the terms first, second, etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, features defined with first, second, etc., may explicitly or implicitly include one or more of that feature. In the description of this invention, unless otherwise stated, multiple means two or more.
[0023] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the terms installation, connection, and linking should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections; they can refer to mechanical connections or electrical connections; they can refer to direct connections or indirect connections through an intermediate medium; and they can refer to the internal communication between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.
[0024] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.
[0025] This invention provides a satellite orbit ECI anchoring rendering method based on osgEarth, comprising the following steps: Step 1, ECI Inertial Frame Static Storage and GMST Rotation Rendering: In the ECI coordinate system, sample the orbital geometry data for half an orbital period before and after the current simulation time, and store it as a point set P_eci. The number of sampling points N is dynamically determined according to the orbital period T_orb. Construct a scene graph node hierarchy structure, which consists of a MapNode (map root node), an independent parent node P_orbit, a MatrixTransform node M_orbit, and orbital line geometry nodes. For each frame, construct a rotation matrix based on the GMST angle calculated at the current simulation time, and map the entire ECI orbital point set to the ECEF world coordinate system. Only update the transformation matrix of the MatrixTransform node M_orbit for each frame. Step 2, Anchor point translation correction: Monitor the deviation between the actual ECI position of the satellite and the orbital geometric center anchor point p_anchor. When the deviation exceeds the recentering threshold, calculate the offset vector and perform an overall vector addition translation on all N points in the point set P_eci. Step 3, Dual-modal orbit display and asynchronous calculation: Supports switching between inertial orbit mode and nadir trajectory mode. Orbit geometry sampling is executed asynchronously in a background thread. Timing coordination between the main thread and the background thread is completed through asynchronous request identifier generation, transmission and matching filtering. Step 4, GPU fragment shader Earth occlusion culling: Turn off depth testing in the rendering pipeline, implement the intersection determination between the viewpoint and the Earth sphere in the custom GLSL fragment shader, calculate the intersection point of the ray from the viewpoint to the current fragment and the Earth sphere, and perform the discard operation when the fragment is behind the Earth. Step 5: Unify the simulation timeline and align it with the spacetime of Earth and Earth: Using the simulation reference time plus the simulation elapsed time as the unified time reference, the satellite SGP4 propagation and the ground entity kinematic model use the same simulation time drive, and all entities are unified to the Earth scene map through the geographic transformation node of osgEarth.
[0026] In a preferred embodiment of the present invention, such as Figure 1 As shown, this method aims to address technical challenges in the rendering of satellite orbit lines in 3D Earth scenes, including orbit line drift, difficulties in handling Earth occlusion by high-orbit satellites, and inconsistencies in spatiotemporal references between Earth and space, in the fields of space situation simulation and digital twins. In these fields, it is necessary to render satellite orbit lines in real-time within a 3D Earth scene and display them collaboratively with ground entities such as vehicles, aircraft, and ships within the same scene. Current technologies typically store orbit sampling points in the ECEF geocentric-ground-fixed coordinate system. Due to Earth's rotation, orbit lines sampled at fixed epochs continuously drift with the Earth's rotation. If SGP4 orbit propagation data is resampled every frame to update the orbit geometry, the CPU overhead of performing tens of thousands of SGP4 propagations per frame for large-scale satellite scenes would make real-time rendering impractical. Furthermore, the orbit lines of high-orbit satellites cross the Earth's surface, and the standard depth buffer lacks sufficient accuracy at long-distance scales, making it difficult to correctly distinguish the relative positions of the orbit lines and the Earth's surface, resulting in penetration artifacts. Simultaneously, the coordinate systems and time references of satellite SGP4 propagation and ground entity kinematic models are inconsistent, lacking a unified simulation time axis driving mechanism.
[0027] To address the aforementioned issues, this method sequentially performs steps such as ECI inertial frame static storage and GMST rotation rendering, anchor point translation correction, dual-modal orbit display and asynchronous computation, GPU fragment shader earth occlusion culling, and unified simulation timeline alignment with spacetime, ultimately achieving drift-free rendering of the osgEarth scene graph and collaborative display of spacetime.
[0028] Step 1: ECI inertial frame static storage and GMST rotation rendering.
[0029] The core idea of this step is to establish a pipeline architecture that stores data once and rotates it every frame, transforming the trajectory drift problem into a fixed matrix update operation. For example... Figures 2 to 3 As shown, the ECI anchoring mechanism consists of a three-layer structure: static ECI orbit geometry storage, MatrixTransform GMST rotation per frame, and anchor point translation.
[0030] In the ECI geocentric inertial coordinate system, the satellite orbit is represented as an ellipse fixed relative to the star, while the osgEarth scene graph is rooted in the ECEF geocentric geofixed coordinate system. To eliminate drift, this step samples the orbital geometry data for half an orbital period before and after the current simulation moment in the ECI coordinate system, storing it as a point set P_eci. This point set contains N sampling points, denoted as P_eci containing p1 to p_N. The number of sampling points N is not a fixed value, but is dynamically determined based on the orbital period T_orb to balance rendering accuracy and computational performance. The formula for calculating the number of sampling points N is as follows: (1); Where T_orb is the orbital period in minutes; K_samp is the sampling density coefficient, representing the number of samples per minute; N_max is the upper limit of the number of sampling points, mainly to prevent excessive sampling points in long-period orbits from causing memory overflow; N_min is the lower limit of the number of sampling points, used to ensure the basic geometric accuracy of short-period orbits. The max function represents taking the maximum value, and the min function represents taking the minimum value. Sampling is performed symmetrically around the current simulation time, and the current ECI position of the satellite is used as the midpoint anchor point of the polyline and forcibly inserted into the middle position of the point set. This anchor point is denoted as p_anchor. This dynamic sampling strategy ensures the geometric accuracy of orbits with different periods and avoids problems caused by oversampling or undersampling.
[0031] Secondly, construct the scene graph node hierarchy. For example... Figure 2As shown in Layer 3, the node hierarchy is as follows: the MapNode (map root node) connects to the independent parent node P_orbit, then to the MatrixTransform node M_orbit, and finally to the track geometry node. The independent parent node P_orbit is mounted to the MapNode as an osg::Group container. Its function is to isolate the large-scale bounding sphere of the track from the calculation of the near and far clipping planes of the main scene camera, preventing the large-scale track from causing the near clipping plane to be too far away and thus culling ground details. The MatrixTransform node M_orbit is a child node of P_orbit, and the track geometry is a child node of M_orbit. This hierarchical structure ensures that the matrix transformation of M_orbit only affects the track and not the main scene.
[0032] After creating a MatrixTransform node M_orbit in the scene graph, a rotation matrix is constructed for each frame based on the GMST angle calculated at the current simulation time. The GMST angle is denoted as... This angle is the Greenwich Mean Sidereal Time, reflecting the Earth's rotation relative to inertial space. The rotation matrix R_z( The calculation formula for ) is as follows: (2); Where cos and sin are the cosine and sine functions, respectively. This matrix represents rotation about the Z-axis. Angle, used to transform the direction of a point in the ECI coordinate system to the ECEF coordinate system.
[0033] Subsequently, the entire ECI orbital point set is mapped to the ECEF world coordinate system. The transformation matrix M is calculated using the following formula: (3); Where T(p_ecef) is the translation matrix of the satellite's current ECEF position, used to translate the orbital geometry from the ECI origin to the satellite's ECEF position. This formula adopts the OSG row vector convention, which actually means first rotating around the Z-axis. The ECI coordinates are converted to ECEF coordinates and then translated to the satellite's current position. Only the M_orbit transformation matrix is updated per frame. The coordinate system mapping is completed using the scene graph's own matrix operation pipeline. Once the orbital geometry data is generated, it resides in the ECI coordinate system and does not change over time. This architecture of storing once and converting per frame means that the main thread only needs to update the rotation matrix once per frame, regardless of the number of satellites, greatly reducing CPU overhead.
[0034] In addition, the ECEF coordinates of the satellite's position are obtained through the following transformation, as shown in the following formula: (4a); (4b); (4c); Where p_eci_x, p_eci_y, and p_eci_z are the x, y, and z coordinate components of the satellite in the ECI coordinate system, respectively; and p_ecef_x, p_ecef_y, and p_ecef_z are the x, y, and z coordinate components of the satellite in the ECEF coordinate system, respectively. Through the above transformation, the spatial relationship between the satellite body and the orbital line remains consistent.
[0035] Step 2: Anchor point translation correction.
[0036] Due to the limitations of floating-point precision and the long-term evolution of the orbit, the statically stored orbital geometric center anchor point p_anchor will deviate from the actual ECI position of the satellite. If not corrected over a long period, the orbital geometry will gradually separate from the satellite's actual position. This step monitors the deviation and corrects it to ensure the geometric accuracy during long-term operation.
[0037] First, monitor the deviation between the satellite's actual ECI position p_eci and the orbital geometric center anchor point p_anchor. Deviation The calculation formula is as follows: (5); The vertical line represents the vector magnitude, which is the Euclidean distance between the two points.
[0038] when Exceeding the re-centering threshold When this occurs, a correction mechanism is triggered. The threshold is reset. Set within the effective precision range of double-precision floating-point, for example Meters. At this point, calculate the offset vector. And perform a global vector addition translation on all N points in the point set P_eci. Offset vector The formula for calculating the translated point p_i' is as follows: ; (6); Where p_i is the i-th sampling point before the translation, and p_i' is the i-th sampling point after the translation. The translation operation involves only N vector additions, and the computational cost is negligible.
[0039] This step employs a dual-path collaborative mechanism. For example... Figure 3 As shown, anchor point translation serves as the fast path, correcting short-term positional deviations; asynchronous full-cycle resampling serves as the slow path, triggered once every T_refresh time, correcting long-term accumulated floating-point errors. During this process, fast path anchor point translation correction is performed; simultaneously, a background thread asynchronously resamples the entire loop every T_refresh time. This dual-path collaboration ensures long-term geometric accuracy and avoids performance degradation caused by frequent resampling.
[0040] Step 3: Dual-modal track display and asynchronous calculation.
[0041] This step supports switching between inertial orbit mode and nadir trajectory mode, and avoids blocking the main thread rendering through asynchronous calculation mechanism.
[0042] like Figure 4 As shown, the dual-mode orbit display includes two forms. In inertial orbit mode, the orbit line is displayed as an ellipse fixed relative to inertial space, rotating as a whole with the Earth's rotation. This mode is suitable for demonstrating the satellite's motion in inertial space. In nadir trajectory mode, geographic coordinates (longitude, latitude, and altitude) are sampled and converted to world coordinates by osgEarth, presenting a figure-eight or sine shape with the Earth's rotation. This mode is suitable for demonstrating the satellite's coverage of the ground.
[0043] like Figure 5 As shown, asynchronous computation coordinates the timing between the main thread and background threads through the generation, transmission, and matching filtering of asynchronous request identifiers. Orbit geometry sampling is executed asynchronously in the background thread, without blocking the main thread's rendering. Each time a mode switch or forced refresh occurs, the main thread increments the asynchronous request identifier `requestId`. When the asynchronous result is returned, the main thread checks whether the request identifier associated with the returned result matches the latest request. If they do not match, the expired result is discarded, ensuring that the interface always displays the latest data consistent with the current selection. This mechanism effectively solves the data race problem in a multi-threaded environment.
[0044] Step 4: GPU fragment shader earth occlusion culling.
[0045] To address the insufficient accuracy of depth buffering caused by high-orbit satellites crossing the Earth's surface, this step performs occlusion culling on the GPU. Geostationary orbit satellites, at an altitude of approximately 35,786 km, have orbits that cross the Earth's surface. Standard depth buffering is insufficiently accurate at long-distance scales, making it difficult to correctly distinguish the relative positions of the orbital line and the Earth's surface, resulting in penetration artifacts. Figure 6 As shown, depth testing is disabled in the rendering pipeline, and the intersection of the line of sight and the Earth sphere is determined in a custom GLSL fragment shader.
[0046] The shader receives the Earth's center c_earth and Earth's radius R_earth as input parameters. c_earth represents the position in the orbital local coordinate system, calculated each frame by the inverse transformation of the MatrixTransform node, ensuring occlusion detection is consistent with the current camera viewpoint. The intersection of the ray from viewpoint o to the current fragment p with the Earth's surface is calculated. The formulas for calculating the ray direction d and the projected length t_c of the line connecting the line of sight to the Earth's center are as follows: ; (7); Here, the dot product symbol represents the vector dot product. t_c represents the projected distance from viewpoint o along the ray direction d to the Earth's center c_earth along the ray direction.
[0047] When t_c > 0, it indicates that the ray direction points towards the Earth's center. In this case, the shortest distance d_min from the line of sight to the Earth's center is calculated. d_min satisfies the following relationship: (8); If d_min 2 <R_earth 2 Then the ray intersects the sphere. The intersection point parameter t_enter is then calculated using the following formula: (9); Here, the radical sign represents the square root operation. t_enter represents the distance from viewpoint o, along the ray direction d, into the Earth's surface.
[0048] When t_enter>0 and t_enter<|po|- At that time, the fragment was located behind Earth and was discarded. To mitigate tangent direction tolerance and prevent misjudgments caused by floating-point precision issues in fragments near the tangent, this step offloads occlusion detection to the GPU fragment shader for parallel execution. In numerous satellite scenes, the overhead of occlusion detection remains constant, only increasing the computational load of the fragment shader's arithmetic logic unit, effectively resolving penetration artifacts caused by insufficient depth buffer precision in high-orbit scenes.
[0049] Step 5: Align the unified simulation timeline with the spacetime of heaven and earth.
[0050] like Figures 7 to 9 As shown, this step adds the simulation elapsed time to the simulation reference time T_ref. To unify the time reference, the satellite SGP4 propagation and ground entity kinematics models use the same simulation time driver, i.e., a unified simulation time axis: T_ref+ The satellite SGP4 propagation and the ground entity kinematic model are connected to the same time reference, which solves the problem of inconsistency between space and time. When the time axis is fast-forwarded / paused, the spatial relative relationship between the satellite and the ground entity remains stable.
[0051] satellite via Perform SGP4 orbital propagation, output ECI coordinates, and rotate them using GMST to obtain ECEF coordinates. Ground entities use... As a time-step driven kinematic model and physics engine, all entities are unified to the Earth scene map through the geographic transformation nodes of osgEarth. This unified timeline mechanism solves the spatiotemporal alignment problem of ground and space entities, ensuring that the alignment of satellite positions with orbital lines remains correct, and that the spatial relative relationships of various types of entities remain stable after performing timeline fast-forward and pause operations.
[0052] In addition, when the camera altitude is below the atmospheric hiding threshold H_hide, the orbital lines are forcibly hidden to avoid visual confusion caused by the orbital lines when observed from within the atmosphere.
[0053] Example 1: Complete Process of ECI Anchoring and GMST Rotation This embodiment is implemented using the osgEarth 3D Earth platform and the libsgp4 open-source orbit propagation library. The parameter settings are as follows: the sampling density coefficient K_samp is... The upper limit of the number of sampling points, N_max, is 360 in inertial orbit mode and 720 in nadir trajectory mode. The lower limit of the number of sampling points, N_min, is 60. The threshold value is then redefined. Pick Meters, slow path resampling interval T_refresh is 5 seconds, atmospheric hiding threshold H_hide is 100km, orbital baseline width is 3.0 pixels, Earth radius R_earth is 6378137m, tangent tolerance Take 50m.
[0054] Step 1: Static storage of ECI coordinates for orbital geometry. Taking a low-Earth orbit satellite with an orbital period of 90 minutes as an example, calculate the number of sampling points N according to formula (1).
[0055] ; In the nadir trajectory mode, N_max is set to 720, and the calculation yields: ; In a background thread, SGP4 propagation is performed symmetrically around the current simulation moment, sampling 360 ECI coordinate points sequentially. The satellite's current ECI position is forcibly inserted as the center anchor point at the 180th position, which is the N / 2th position. The point set P_eci is stored in memory and will not be resampled subsequently.
[0056] Step 2: Scene Graph Node Construction. Create nodes in the osgEarth scene graph according to the hierarchical structure described above. First, create an independent parent node P_orbit and attach it to MapNode to isolate the large-scale bounding sphere of the orbital lines from the near and far clipping plane calculations of the main scene camera. Then, create a MatrixTransform node M_orbit under P_orbit. The orbital geometry is a child node of M_orbit. The node hierarchy is as follows: the MapNode root node connects to the independent parent node P_orbit, connects to the M_orbit transform node, and connects to the orbital geometry node. Ensure that the matrix transformation of M_orbit only affects the orbital lines and not the main scene.
[0057] Step 3: GMST rotation update per frame. The main thread calculates the GMST angle at the current simulation moment in each frame. Let the current situation be... -45 degrees corresponds to 3 hours before Greenwich Mean Time. Substitute into formula (2) to construct the rotation matrix. cos(-45 degrees) = 0.7071, sin(-45 degrees) = -0.7071. Substitute into formula (3), using the satellite's current ECEF position as the translation component, to obtain the transformation matrix of M_orbit. Set M as the matrix of the M_orbit node. The main thread only performs one trigonometric function calculation and one matrix multiplication, which is independent of the number of satellites in the scene. When switching back to the ground reference frame for observation, the orbital line rotates as a whole with the Earth's rotation, but the relative position of the orbital line and the satellite body remains unchanged.
[0058] Step 4: Anchor point translation triggering and execution. The system continuously monitors the deviation between the satellite's actual ECI position and the center anchor point p_anchor. After running for several minutes, the satellite moves along its orbit to adjust the deviation. achieve Meters. The anchor point p_anchor is the position at the sampling time. After the satellite continues to move, the ECI position naturally deviates from the anchor point. At this time, the deviation exceeds the recentering threshold. Meters. Substitute into formula (5) to calculate the deviation, and substitute into formula (6) to calculate the offset vector. Vector addition and translation are performed on all 360 points in P_eci. The total computational cost of translation is 360 three-dimensional vector additions, which takes microseconds. When the satellite continues to operate until the T_refresh=5-second interval, the background thread triggers asynchronous full-cycle resampling, resampling 360 ECI points centered on the current simulation time to correct long-term accumulated errors.
[0059] Step 5: Dual-mode switching and asynchronous calculation verification. The dual-mode orbit display from Step 3 switches between inertial orbit mode and nadir trajectory mode. In inertial orbit mode, the orbital line is a fixed ellipse; in nadir trajectory mode, the orbital path is sampled using geographic coordinates and presents a figure-eight shape as the Earth rotates.
[0060] The asynchronous computation process is as follows: Each time track resampling is triggered, the main thread increments the global request identifier `requestId`, which is a monotonically increasing 64-bit unsigned integer. The `requestId`, along with the current track display mode type, TLE data, current simulation time, and other parameters, is packaged and passed to the background thread. After the background thread completes SGP4 propagation sampling, it returns the `requestId`, track display mode type, and computation result. The main thread performs two checks sequentially when consuming the result. The first check checks if the returned `requestId` is equal to the latest request identifier `m_orbitRequestId`. If they are not equal, it means an update request has already been issued, and the expired result is discarded. The second check checks if the returned track display mode type is consistent with the current display mode. If they are inconsistent, it means the user switched modes during asynchronous computation, and this result is also discarded. After both checks pass, the computation result is applied to the track geometry nodes in the scene graph.
[0061] When rapidly switching between two modes, for example, switching five times consecutively with 200ms intervals, each switch increments the requestId and sends a new asynchronous request. Results of previously incomplete requests are automatically discarded due to mismatched requestIds. The interface always displays the latest data consistent with the current selection, without any old asynchronous results overwriting new results or mode misalignment.
[0062] Example 2: Complete Shading Process for GPU Occlusion Culling This embodiment covers steps four and five of the invention.
[0063] Step 1: Shader Input Parameter Calculation. Taking a geostationary satellite at an altitude of 35786 km as an example, the main thread calculates the Earth's center position (c_earth) in the local orbital coordinate system and the Earth's radius (R_earth = 6378137 m, i.e., the WGS-84 equatorial radius) each frame using the inverse transformation of the MatrixTransform node. These two values are then passed as a uniform to the shader. The viewpoint (o) is obtained from the position of the osgEarth camera in the world coordinate system through the same inverse transformation.
[0064] Step 2: Occlusion determination in the fragment shader. The camera is located on the back side of the Earth to observe the GEO satellite orbit line. For each fragment p on the orbit line, substitute it into formula (7) to calculate t_c. Let |c_earth-o|=42000000m, which is the distance from the viewpoint to the Earth's center. R_earth=6378137m, then R_earth 2 =4.068×10 13 m 2 .
[0065] For fragments that are not occluded in front of the Earth, d points in the direction away from the Earth, t_c<0, the shader skips the occlusion determination and outputs the fragment color normally.
[0066] For occluded fragments behind the Earth, let t_c=40000000m, substitute into formula (8) to calculate d_min 2 .
[0067] ; Since d_min 2 =1.64×10 14 >4.068×10 13 , that is d_min>R_earth, the ray grazes past the Earth and is not determined as occluded.
[0068] In another case, let t_c=41900000m.
[0069] ; Since 8.39×10 12 <4.068×10 13 , that is d_min<R_earth, the ray passes through the Earth sphere. Substitute into formula (9) to calculate t_enter.
[0070] ; The distance of the current fragment |p-o|=4.21×10 7 m>t_enter, execute discard to discard this fragment.
[0071] Step 3: Comparative verification. When the custom shader is turned off and falls back to the standard depth test in the same scene, the depth buffer has insufficient precision at the long-distance scale, and the orbit line penetrates the Earth's surface to produce artifacts. After enabling the custom shader, the penetration artifacts are eliminated. The line segments of the orbit line on the back side of the Earth are correctly discarded, and the orbit line segments near the horizon present a smooth fade-out effect. Depth testing and lighting are turned off in the rendering pipeline settings, and visibility is completely determined by the geometric judgment in the fragment shader.
[0072] Step 4: Verify the alignment of spacetime and time. Step 5 uses the unified simulation timeline, which is the simulation reference time T_ref plus the simulation elapsed time. Using a time reference, the propagation of the SGP4 satellite is synchronized with the kinematic model of ground entities. When satellites, ground vehicles, and aircraft are deployed simultaneously in the same scene, and a unified simulation timeline is set, the alignment of satellite positions with orbital lines remains correct. After performing timeline fast-forward and pause operations, the spatial relative relationships of various types of entities remain stable. When the camera altitude drops below H_hide=100km, the orbital lines are forcibly hidden; when they rise above 100km, they are redisplayed to avoid visual confusion caused by orbital lines when observing from within the atmosphere.
[0073] Compared with existing technologies, the satellite orbit ECI anchoring rendering method based on osgEarth provided by this invention has the following significant advantages: First, this invention completely solves the problems of orbital drift and rendering performance bottlenecks in large-scale satellite scenes by establishing a pipeline architecture that combines static storage of orbital geometry in the ECI inertial frame with GMST rotation of the MatrixTransform node per frame—storing once and rotating every frame. Existing technologies typically store orbital points in the ECEF coordinate system, leading to drift with Earth's rotation, or re-sampling SGP4 every frame, resulting in excessive CPU computational overhead. This invention performs orbital geometry sampling once and statically stores it in the ECI coordinate system. Utilizing the scene graph's own matrix operation pipeline, dynamic coordinate system mapping is achieved by performing a rotation matrix update only once per frame. This architecture reduces the computational complexity of orbital rendering from a linear relationship proportional to the number of satellites to a constant level, greatly freeing up CPU resources and enabling smooth real-time rendering frame rates even in large-scale scenes containing hundreds of satellites. Simultaneously, the number of sampling points is dynamically determined by a formula, balancing geometric accuracy and memory usage for different orbital periods.
[0074] Secondly, this invention innovatively designs an anchor point translation correction mechanism, effectively balancing real-time rendering performance with long-term geometric accuracy. Due to the limitations of computer floating-point precision, long-running simulations can lead to accumulated errors between the statically stored orbital geometric center and the actual satellite position. This invention employs a fast-slow dual-path collaborative strategy: within a small timescale, fast path translation correction is performed by monitoring the deviation and executing batch vector addition, an operation with extremely low computational overhead; within a long timescale, slow path correction is performed by asynchronously triggering full-cycle resampling via a background thread. This mechanism avoids the performance loss caused by frequent SGP4 propagation and fundamentally eliminates long-term accumulated floating-point errors, ensuring the absolute accuracy of orbital geometry in long-term simulation tasks.
[0075] Third, the asynchronous computation and requestId matching filtering mechanism proposed in this invention significantly improves the system's response speed and data consistency. During dual-modal track display switching or forced refresh, traditional synchronous computation can block the main thread, causing interface lag. This invention offloads the time-consuming SGP4 sampling task to a background thread for asynchronous execution and introduces a monotonically increasing requestId identifier. When the user frequently switches modes during asynchronous computation, the main thread can accurately identify and automatically discard expired computation results, ensuring that the interface always renders the latest data consistent with the current user operation. This not only eliminates the risk of data contention in a multi-threaded environment but also guarantees extremely smooth user interaction.
[0076] Fourth, this invention proposes an Earth occlusion culling scheme based on GPU fragment shaders, perfectly solving the problem of insufficient depth buffer accuracy in high-orbit satellite rendering. In large-scale scenes at long distances, such as geostationary orbit, traditional Z-buffer depth testing is prone to Z-fighting artifacts, where orbital lines penetrate the Earth's surface, due to accuracy loss. This invention disables depth testing in the rendering pipeline and writes the occlusion determination logic into a custom GLSL fragment shader in the form of analytical geometric formulas. By calculating the precise intersection parameters between the line of sight and the Earth's sphere, fragment-level occlusion culling is performed in parallel on the GPU. This scheme not only completely eliminates penetration artifacts and achieves a smooth fading effect near the horizon, but also evenly distributes the computational burden across the GPU's arithmetic logic units, keeping the occlusion determination overhead constant and resulting in extremely high rendering efficiency.
[0077] Finally, this invention achieves high-precision spatiotemporal alignment of ground and space entities by constructing a unified simulation timeline. Addressing the pain point of inconsistencies between satellite orbit propagation and the coordinate systems and time bases of ground entity kinematic models, this invention uses the sum of the simulation reference time and the elapsed time as the sole driving source, uniformly scheduling SGP4 propagation and physics engine stepping. By combining the hierarchical combination of MatrixTransform and geographic transformation nodes in the scene graph, it ensures that the spatial relative relationships of the satellite body, orbit lines, and ground vehicles and spacecraft remain strictly consistent throughout the fast-forwarding, pause, or normal elapsed timelines. Furthermore, by incorporating a dynamic camera altitude hiding mechanism, it effectively avoids visual clutter caused by orbit lines from low-altitude perspectives, significantly enhancing the immersiveness and professionalism of the 3D Earth scene.
[0078] In summary, this invention has undergone in-depth optimization in multiple dimensions, including rendering architecture, precision control, concurrent processing, graphics algorithms, and system collaboration, providing an efficient, accurate, and stable satellite orbit visualization solution for the fields of aerospace situation simulation and digital twins.
[0079] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.
Claims
1. A satellite orbit ECI anchoring rendering method based on osgEarth, characterized in that, Includes the following steps: Step 1, ECI Inertial Frame Static Storage and GMST Rotation Rendering: In the ECI coordinate system, sample the orbital geometry data for half an orbital period before and after the current simulation moment at once, and store it as a point set P_eci. The number of sampling points N is dynamically determined according to the orbital period T_orb. Construct a scene graph node hierarchy structure, which consists of MapNode (map root node), independent parent node P_orbit, MatrixTransform (transformation node M_orbit), and orbital line geometry nodes. Each frame constructs a rotation matrix based on the GMST angle calculated at the current simulation time, mapping the entire ECI orbit point set to the ECEF world coordinate system. Each frame only updates the transformation matrix of the MatrixTransform transformation node M_orbit. Step 2, Anchor point translation correction: Monitor the deviation between the actual ECI position of the satellite and the orbital geometric center anchor point p_anchor. When the deviation exceeds the recentering threshold, calculate the offset vector and perform an overall vector addition translation on all N points in the point set P_eci. Step 3, Dual-modal orbit display and asynchronous calculation: Supports switching between inertial orbit mode and nadir trajectory mode. Orbit geometry sampling is executed asynchronously in a background thread. Timing coordination between the main thread and the background thread is completed through asynchronous request identifier generation, transmission and matching filtering. Step 4, GPU fragment shader Earth occlusion culling: Turn off depth testing in the rendering pipeline, implement the intersection determination between the viewpoint and the Earth sphere in the custom GLSL fragment shader, calculate the intersection point of the ray from the viewpoint to the current fragment and the Earth sphere, and perform the discard operation when the fragment is behind the Earth. Step 5: Unify the simulation timeline and align it with the spacetime of Earth and Earth: Using the simulation reference time plus the simulation elapsed time as the unified time reference, the satellite SGP4 propagation and the ground entity kinematic model use the same simulation time drive, and all entities are unified to the Earth scene map through the geographic transformation node of osgEarth.
2. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step one, the number of sampling points N is dynamically determined based on the orbital period T_orb using the following formula: ; Wherein, K_samp is the sampling density coefficient, N_max is the upper limit of the number of sampling points, and N_min is the lower limit of the number of sampling points; sampling is performed symmetrically around the current simulation time, and the current ECI position of the satellite is used as the midpoint anchor point of the polyline and forcibly inserted into the middle position of the point set.
3. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step one, a rotation matrix is constructed for each frame based on the GMST angle calculated at the current simulation time, mapping the entire ECI orbit point set to the ECEF world coordinate system. Specifically, this includes: Construct the rotation matrix R_z(θ_gmst) based on the GMST angle θ_gmst calculated at the current simulation time: ; The transformation matrix M for mapping the entire ECI orbit point set to the ECEF world coordinate system is calculated as follows: ; Where T(p_ecef) is the translation matrix of the satellite's current ECEF position.
4. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step two, the deviation between the satellite's actual ECI position and the orbital geometric center anchor point p_anchor is monitored. When the deviation exceeds the recentering threshold, the offset vector is calculated, and a global vector addition translation is performed on all N points in the point set P_eci. Specifically, this includes: Calculate the deviation: ; When δ exceeds the recentering threshold ε, calculate the offset vector: ; Perform global vector addition and translation on all N points in the point set P_eci: ; Where p_i is the i-th sampling point before translation, and p_i' is the i-th sampling point after translation.
5. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step three, the timing coordination between the main thread and the background thread is completed through asynchronous request identifier generation, transmission, and matching filtering. Specifically, this includes: Each time a mode switch or forced refresh occurs, the main thread increments the asynchronous request identifier requestId and passes the requestId along with the current track display mode type to the background thread; After the background thread completes the sampling, it returns the requestId, track display mode type, and calculation result together. When consuming results, the main thread performs two checks in sequence: the first check checks whether the returned requestId is equal to the latest request identifier, and the second check checks whether the returned track display mode type is consistent with the current display mode; after both checks pass, the calculation results are applied to the track geometry nodes of the scene graph.
6. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step four, the intersection determination between the viewpoint and the Earth sphere is implemented in the custom GLSL fragment shader, specifically including: The shader takes the Earth's center c_earth and Earth's radius R_earth as input parameters and calculates the ray direction: ; And the projected length of the line of sight to the center of the Earth: ; When t_c>0, the shortest distance d_min from the line of sight to the center of the Earth satisfies: ; If d_min 2 <R_earth 2 Then calculate the intersection parameters of the points entering the sphere: ; When t_enter > 0 and t_enter < |po| - ε_t, a discard operation is performed; where ε_t is the tangent direction tolerance.
7. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step five, the simulation reference time plus the simulation elapsed time is used as a unified time reference. The satellite SGP4 propagation and ground entity kinematic models use the same simulation time drive, specifically: Satellite passes by: ; SGP4 orbital propagation is performed, and the output ECI coordinates are rotated by GMST to obtain ECEF coordinates; the ground entity uses Δt_sim as the time step to drive the kinematic model and physics engine. In addition, when the camera altitude is below the atmospheric hiding threshold H_hide, the orbital lines are forcibly hidden.
8. The satellite orbit ECI anchoring rendering method based on osgEarth according to claim 1, characterized in that, In step three, the dual-mode orbital display includes: In inertial orbit mode, the orbit line appears as an ellipse fixed relative to inertial space, rotating as a whole with the Earth's rotation. In the nadir trajectory mode, the sampled geographic coordinates are converted from osgEarth to world coordinates and appear in a figure-eight or sine shape as the Earth rotates.
Citation Information
Patent Citations
Satellite multi-tense visualization method based on SGP4 orbit model and database
CN120804206A
Satellite coordinate adaptive conversion method and system
CN121877022A