Aircraft physical simulation and attitude control method based on osgEarth

Abstract

The application discloses a method for aircraft physical simulation and attitude control based on osgEarth and relates to the technical field of computer simulation.The method comprises the following steps: fixed-wing dynamics is reduced to horizontal and vertical velocity independent calculation, and a segmented lift model is adopted;target attitude angles are calculated by a first-order exponential smoothing and a coordinated turning mechanism in the kinematics layer;target attitude angles are constructed into a quaternion in a specified order and written into a physical engine in a one-way manner, and rigid body angular velocity is frozen;after a physical step, the current attitude angle saved by the kinematics layer is directly taken to update an OSG scene graph, and a quaternion readback is skipped.The application changes the attitude update path from two-way conversion to one-way writing, eliminates the precision loss of inverse trigonometric functions and the gimbal lock problem, avoids error accumulation drift, significantly reduces the calculation complexity, and meets the performance and precision requirements of large-scale real-time situation deduction.

Description

Technical Field

[0001] This invention belongs to the field of computer simulation and digital twin technology, and in particular relates to a method for physical simulation and attitude control of aircraft based on osgEarth. Background Technology

[0002] In digital twin potential simulation systems, it is often necessary to introduce a general-purpose physics engine into a 3D digital earth scene (such as the osgEarth platform) to perform highly realistic aircraft simulations. However, different engines differ in how they store and represent attitude data: general-purpose physics engines (such as Bullet) typically use quaternions to store rigid body attitudes to avoid singularity issues; while 3D scene graph engines (such as OpenSceneGraph, or OSG) typically use heading-pitch-roll (HPR, or Euler angles) to represent node attitudes.

[0003] To achieve synchronization between physical simulation and scene rendering, existing solutions typically perform a bidirectional conversion process every frame, from Euler angles to quaternions and from quaternions back to Euler angles. Specifically, the system converts the kinematic target (HPR) into quaternions and writes them into the physics engine. After the physics engine performs dynamic integration, it reads back the quaternions and converts them back to HPR to update the scene graph.

[0004] However, the aforementioned existing technical solutions have the following significant defects and shortcomings in practical applications: First, this bidirectional conversion process suffers from severe precision loss and numerical instability. The conversion from quaternions to HPR involves the extraction of rotation matrix elements and inverse trigonometric function operations such as arcsine (asin) and arctangent (atan2). Taking the ZXY rotation order as an example, the pitch angle is extracted using the arcsine function, while the heading and roll angles are extracted using the arctangent function. When the pitch angle approaches ±90°, the cosine value approaches zero, causing the denominator of the arctangent function to approach zero, making the calculation of the heading and roll angles highly prone to divergence. Furthermore, the rounding error of double-precision floating-point numbers is significantly amplified in the inverse trigonometric function operations.

[0005] Secondly, existing technologies struggle to overcome the gimbal lock problem. When the pitch angle reaches ±90°, the yaw and roll axes coincide in space, causing the system to lose a rotational degree of freedom, or HPR (uniqueness). At the edge of gimbal lock, tiny floating-point rounding errors, amplified by derivatives, can produce invalid or outlier values, severely impacting the reliability of the simulation.

[0006] Third, the aforementioned conversion errors accumulate frame by frame during continuous roll and other maneuver control, resulting in a non-negligible drift in the heading angle, which cannot meet the requirements of high-precision situation simulation.

[0007] Finally, existing flight simulation frameworks typically employ complete six-degree-of-freedom rigid body dynamics equations and obtain parameters through aerodynamic coefficient tables. This approach has extremely high computational complexity, and while suitable for professional flight training scenarios, it struggles to meet the stringent computational performance requirements of large-scale, real-time situation simulation systems. Furthermore, current technologies lack mechanisms for directly controlling attitude through kinematic layers within a specific physics engine integration environment, and also lack dedicated simplified aerodynamic models for situation simulation scenarios.

[0008] Therefore, there is an urgent need to develop an aircraft physics simulation and attitude control method suitable for the osgEarth platform and Bullet physics engine to solve the problems of low attitude synchronization accuracy, susceptibility to gimbal lock, and excessive computational overhead. Summary of the Invention

[0009] In view of this, the present invention aims to propose an aircraft physics simulation and attitude control method based on osgEarth to solve the technical defects of existing digital twin potential inference systems, such as the accuracy loss caused by bidirectional conversion of HPR (Euler angles) and quaternions, gimbal lock risk and frame-by-frame error accumulation drift, and the high computational complexity of the complete six-degree-of-freedom rigid body dynamics model, which makes it difficult to meet the performance requirements of large-scale real-time situation inference.

[0010] This invention provides a method for aircraft physics simulation and attitude control based on osgEarth, including the following steps: Step 1, velocity separation and simplified aerodynamic model calculation: reduce the six-degree-of-freedom rigid body dynamics of the fixed wing to two independent one-dimensional problems of horizontal velocity and vertical velocity, and calculate the combined velocity, lift, thrust, drag and new velocity; Step 2, kinematic attitude calculation: The target attitude angle is calculated using first-order exponential smoothing, and the yaw input is derived through a coordinated turning mechanism; Step 3, One-way writing of kinematics layer and freezing of angular velocity: The target attitude angles calculated by the kinematics layer are constructed as quaternions according to the specified rotation order, written into the physics engine and the rigid body angular velocity is frozen. Step 4, Synchronous Update Layer Skips Quaternion Readback: After the physics step is completed, the current heading angle, pitch angle and roll angle saved in the kinematics layer are directly retrieved to update the OSG scene map, without going through the quaternion to Euler angle readback inverse operation.

[0011] Furthermore, in step one, the velocity vector is synthesized. The calculation formula is: ; In the formula, The magnitude of the horizontal velocity. Vertical velocity, This is the heading angle.

[0012] Furthermore, in step one, a segmented lift model is adopted, with the stall speed... As a dividing line: When horizontal speed At that time, the base lift ; when At that time, the base lift ; Additional lift ,in ; Total lift ; In the formula, For the mass of the aircraft, It is the acceleration due to gravity. The pitch angle, This is the velocity factor.

[0013] Furthermore, in step one, the lift force is decomposed into a vertical component. and horizontal backward component ; Thrust decomposed into and ; resistance ; New speed calculation is , ; In the formula, For total lift, The pitch angle, For vertical components, The horizontal backward component, For maximum thrust, This is the throttle coefficient. For horizontal thrust, For vertical thrust, The drag coefficient, Calculate speed for resistance. As a resistance multiplier, As resistance, For a new level of speed, For horizontal acceleration, For time step, For the new vertical velocity, This is the vertical acceleration.

[0014] Furthermore, in step two, the formula for the smooth transition of the target attitude angle is: ; The formula for coordinated turn yaw input is: ; In the formula, For the updated attitude angles, The attitude angle of the previous frame. For the target attitude angle, To approximate the rate coefficient, For time step, For yaw angle increment, The current roll angle, For the maximum roll angle, To coordinate the turning gain coefficient, This is the maximum yaw rate.

[0015] Furthermore, in step three, the target attitude angle Quaternions are constructed by rotating from Z to X to Y. The order of quaternion multiplication is as follows: ; The components are constructed as follows: , , ; Freeze the rigid body angular velocity to ; In the formula, For the target quaternion, For the heading quaternion, For the quaternion of pitch, For rolling quaternions, For heading angle, The pitch angle, For roll angle, , , These are the radian values ​​corresponding to the heading angle, pitch angle, and roll angle, respectively. ω is the rigid body angular velocity.

[0016] Furthermore, in step four, skipping the quaternion backread specifically means: not executing the inverse formula from quaternions to Euler angles. , , ; In the formula, The pitch angle, For heading angle, For roll angle, For rotation matrix, , , , , Rotation matrices The elements in the corresponding row and column.

[0017] Based on the same concept, the present invention also provides a computer-readable storage medium having a computer program stored thereon, characterized in that the computer program, when executed by a processor, implements an osgEarth-based method for aircraft physical simulation and attitude control.

[0018] Compared with existing technologies, the osgEarth-based aircraft physics simulation and attitude control method of the present invention has the following advantages: (1) The accuracy loss and gimbal lock risk caused by attitude readback are completely eliminated. This invention reduces the attitude update path from two bidirectional conversions in the traditional scheme to one unidirectional conversion by directly calculating the attitude angle at the kinematic layer, writing and freezing the angular velocity in one direction at the physical layer, and skipping the quaternion readback at the synchronization layer. This method avoids the accuracy loss of arcsine and arctangent operations near the gimbal lock in the quaternion to Euler angle readback. The attitude accuracy depends only on the unidirectional quaternion construction with double-precision floating-point precision, fundamentally solving the problems of numerical instability and cumulative drift of heading angle, and significantly improving the simulation accuracy of situation inference.

[0019] (2) Significantly reduced computational complexity, meeting the requirements of large-scale real-time simulation. This invention reduces the six-degree-of-freedom rigid body dynamics to two independent one-dimensional problems: horizontal velocity and vertical velocity. The lift calculation is driven by the stall velocity as the boundary through a piecewise lift model. The lift calculation only involves algebraic operations without complex iterative integration. The computational load per frame is significantly reduced compared to the complete six-degree-of-freedom model, effectively meeting the stringent requirements of large-scale, real-time situation simulation systems for computational performance.

[0020] (3) Achieving smooth attitude transition without overshoot and highly stable attitude control. This invention achieves attitude transition without overshoot at the kinematic level through a first-order exponential smoothing algorithm and automatically derives the yaw input through a coordinated turning mechanism. At the same time, the physical layer does not drive angular velocity integration through torque, and the attitude control is completely dominated by the kinematic layer, so that the simulation process is not affected by the integral error of the physical engine, thus improving the realism and controllability of the aircraft maneuver simulation. Attached Figure Description

[0021] 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 A schematic diagram of the osgEarth-based physical simulation and attitude control method for aircraft provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of velocity separation and synthesis provided for an embodiment of the present invention, wherein the left side is a horizontal top view and the right side is a vertical cross-sectional side view; Figure 3This is a schematic diagram comparing the attitude update paths of the conventional solution and the solution provided in this embodiment of the invention; Figure 4 This is a schematic diagram of the segmented lift model curve provided in an embodiment of the present invention; Figure 5 This is a schematic diagram of the coordinate system rotation order and quaternion construction provided in an embodiment of the present invention; Figure 6 This is a schematic diagram illustrating the gimbal deadlock principle and the inverse solution principle from quaternions to Euler angles provided in an embodiment of the present invention. Detailed Implementation

[0022] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0023] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, 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, a feature 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, "a plurality of" means two or more.

[0024] 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 a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art will understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0025] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0026] This invention provides a method for physical simulation and attitude control of aircraft based on osgEarth, comprising the following steps: Step 1, velocity separation and simplified aerodynamic model calculation: reduce the six-degree-of-freedom rigid body dynamics of the fixed wing to two independent one-dimensional problems of horizontal velocity and vertical velocity, and calculate the combined velocity, lift, thrust, drag and new velocity; Step 2, kinematic attitude calculation: The target attitude angle is calculated using first-order exponential smoothing, and the yaw input is derived through a coordinated turning mechanism; Step 3, One-way writing of kinematics layer and freezing of angular velocity: The target attitude angles calculated by the kinematics layer are constructed as quaternions according to the specified rotation order, written into the physics engine and the rigid body angular velocity is frozen. Step 4, Synchronous Update Layer Skips Quaternion Readback: After the physics step is completed, the current heading angle, pitch angle and roll angle saved in the kinematics layer are directly retrieved to update the OSG scene map, without going through the quaternion to Euler angle readback inverse operation.

[0027] This method is applied to a digital twin potential extrapolation system, aiming to solve the attitude synchronization accuracy and computational complexity problems encountered when introducing a general physics engine into a 3D digital earth scene for aircraft simulation. In this embodiment of the invention, osgEarth serves as the 3D digital earth scene engine, using heading, pitch, and roll (HPR Euler angles) to represent node attitudes; Bullet serves as the general physics engine, using quaternions to store rigid body attitudes.

[0028] In this embodiment, as Figure 1 As shown, the overall process of this method sequentially executes velocity separation and simplified aerodynamic model calculation, kinematic layer attitude calculation, kinematic layer unidirectional writing and angular velocity freezing, and synchronous update skipping quaternion readback. The specific steps are described in detail below with reference to the accompanying figures.

[0029] Step 1: Velocity separation and simplified aerodynamic model calculation like Figure 2 As shown, this step reduces the six-DOF rigid body dynamics of the fixed wing to horizontal velocity. and vertical velocity Two independent one-dimensional problems. The initial state of the aircraft includes... , , , , The initial value, where For heading angle, The pitch angle, This is the roll angle. The horizontal velocity along the heading angle... Direction is driven solely by thrust and drag, without decomposing the gravitational component. Vertical velocity accumulates independently, driven by lift and gravity. This dimensionality reduction significantly lowers computational complexity, meeting the performance requirements of large-scale real-time situational simulations.

[0030] Composite velocity vector Calculate using the following formula: ; In the formula, The magnitude of the horizontal velocity. For vertical velocity, the ENU coordinate system is used here, which is the East, North, and Sky coordinate system, where the Z-axis is positive upwards; This is the heading angle, with 0° representing north and increasing clockwise. For example... Figure 2 The top view of the horizontal plane on the left side shows the horizontal velocity. along the heading angle Direction decomposed into eastward component and northward component .like Figure 2 The resultant velocity is shown in the vertical cross-sectional side view on the right side of the middle section. From horizontal components With vertical component Orthogonal synthesis.

[0031] like Figure 4 As shown, the segmented lift model uses the stall velocity... To define the boundaries, lift calculations are divided into two sections: the normal flight zone and the stall zone. Basic lift. The calculation is as follows: ; ; In the formula, For the mass of the aircraft, The acceleration due to gravity is taken as 9.8 m / s². 2 When the horizontal speed is not lower than the stall speed, the basic lift equals the weight, which is within the normal flight zone; when it is lower than the stall speed, the lift decreases by the square of the speed ratio, simulating the sharp drop in lift caused by airflow separation above the wing, which is within the stall zone.

[0032] Additional lift It is proportional to the velocity factor and the sine of the pitch angle, and is calculated as follows: ; In the formula, The velocity factor is squared to simulate the square effect of dynamic pressure, with an upper limit of 1.5 to prevent lift divergence during overspeed. This is the pitch angle, with "head up" being positive. Total lift. Total lift with respect to pitch angle and velocity factor The explicit function is .

[0033] Under normal flight conditions ( That is, the pitch angle is greater than the negative critical angle, usually At an angle of 5° to 10°, the lift direction is upward ( , ), clamp To eliminate numerical noise caused by small disturbances (typical disturbance magnitude is about 10). -14 N). When the aircraft performs a negative angle-of-attack maneuver (push-up dive), ) or flying backwards ( )hour, Follow The sign is reversed, resulting in a significantly negative value: in the example parameters ( , )Down, and When s = 1.5, , ,at this time This negative lift is a mathematical expression of the physical fact that the wing generates downward aerodynamic force during a pushrod dive. No force is applied under this condition. Clamping, lift can take negative values. The selection of 5° to 10° is based on the fact that this range lies between the upper limit of the normal descent pitch angle (-5°) and the lower limit of the typical pitch angle for a push-button dive (-15°), providing ample transition space for switching from normal flight to negative lift maneuvers. The specific implementation follows the above... The function's conditional statement enables / disables dynamic enable / disables clamping conditions.

[0034] When lift is decomposed in the ENU coordinate system (East-North-Sky, X / Y / Z), it is determined by the pitch angle. and roll angle Together, they determine its spatial orientation. The direction of lift L is along the Z_B axis of the body coordinate system (perpendicular to the wing plane and upward), and its three-dimensional projection in the world ENU coordinate system is obtained by multiplying the rotation matrices of the Z→X→Y internal rotation: ,in , , These are around the Z-axis (heading) ), around the X' axis (pitch angle) ), around the Y'' axis (roll angle) The rotation matrix of ).

[0035] In the absence of rolling ( In the symmetrical flight of ), the projection degenerates into The lift force is decomposed into a vertical component in the ENU coordinate system. and the horizontal backward component in the opposite direction of the world coordinate system heading : ; ; In the formula The direction is fixed in the opposite direction of the heading in the world ENU coordinate system (-heading direction vector), that is... ,in This refers to the heading angle. When the aircraft has a roll angle... (Roll) When the lift force is applied, it produces an additional lateral component in the horizontal plane, the magnitude of which is obtained by expanding the X component of the complete projection matrix. The direction is perpendicular to the heading direction (in the horizontal plane, along the right side of the heading direction, orthogonal to the heading direction vector in the horizontal plane). This lateral component provides the centripetal acceleration. In steady-state coordinated turning ( , (The restoring process is to the standard coordinated turning formula) In the coordinated turn of the embodiment, the lateral component is given by formula (10) in a kinematic equivalent manner, ensuring the physical consistency of the turning mechanics.

[0036] Thrust from maximum thrust and throttle coefficient Decide, In the ENU coordinate system, it is decomposed into: ; ; drag coefficient Deduced from the maximum speed. At the maximum speed During horizontal flight, thrust and drag are balanced, that is... ,have to The resistance formula is: ; In the formula, As a drag multiplier, drag increases during nose-up or nose-down maneuvers. Maximum pitch angle; When the ground stage, i.e., the relative height to the ground surface, is less than 5 meters, take... Horizontal net force acceleration Limited to Within the interval, This represents the maximum horizontal acceleration. Vertical acceleration... In the formula, vertical damping It only works when rising. , is the vertical damping ratio.

[0037] The new speed is calculated as follows: ; ; The vertical velocity constraint distinguishes between the following two independent conditions: assuming the aircraft's altitude relative to the ground surface... (a) Ground effect judgment conditions ( This is only used to switch the drag calculation model; within this range, air drag is mainly composed of the horizontal component. Pick Rather than synthesis speed It does not involve the (a) Any clamping. (b) Ground contact constraint conditions And it is only applied when the aircraft's wheel load height is zero (actually in contact with the ground). To prevent aircraft from physically penetrating the ground. Throughout the entire height range, Unconstrained by any altitude greater than zero, the aircraft can descend freely according to aerodynamic equations to complete the landing.

[0038] Step 2: Kinematic Attitude Calculation The smooth transition of the target attitude angle adopts first-order exponential smoothing: ; In the formula, For roll angle or pitch angle, For the updated attitude angles, The attitude angle of the previous frame. For the target attitude angle, To approximate the rate coefficient, is the time step. This formula ensures that attitude changes are free from overshoot when there are abrupt changes in the control input.

[0039] The coordinated yaw input is driven by the roll angle: ; In the formula, For yaw angle increment, The current roll angle, For the maximum roll angle, To coordinate the turning gain coefficient, This is the maximum yaw rate.

[0040] Step 3: Kinematic layer dominant attitude, i.e., unidirectional writing and angular velocity freezing. like Figure 5 As shown, the target attitude angles calculated by the kinematic layer are... Quaternions are constructed using the rotation sequence Z to X to Y. This sequence refers to first rotating the yaw angle around the world Z-axis, then rotating the pitch angle around the rotated X'-axis, and finally rotating the roll angle around the rotated Y''-axis. The quaternion multiplication order is as follows: Each component is constructed using axis angles: ; ; ; In the formula, For the heading quaternion, For the quaternion of pitch, For rolling quaternions, For the target quaternion. This follows the OSG coordinate system's negative sign convention; ; It should be noted that the OSG coordinate system is defined as X for right, Y for forward, Z for up, with 0° heading as the positive Y-axis direction and clockwise as positive.

[0041] Will Set the rotational component of the rigid body world transformation matrix and immediately freeze the rigid body angular velocity. The physics engine no longer generates attitude changes through angular velocity integration; attitude control is now entirely delegated to the kinematic layer.

[0042] Angular velocity freeze is persisted through the rigid body angular velocity factor interface of the Bullet physics engine: call angular velocity factor of rigid body Set as a zero vector. Bullet's discrete dynamics world stepSimulation processes the rigid body states in the motion integration phase (integrateTransforms step) of each frame in the following order: (Component-wise product); Linear velocity integral: ; Angular velocity to quaternion increment conversion dq; Attitude update: .when At this time, the effective angular velocities along each axis are always zero vectors, and the quaternion increment dq is always a unit quaternion (1, 0, 0, 0). Therefore... The pose is maintained in all subsequent frames according to the kinematics layer writing. Unchanged. The linear velocity integral path does not pass through... The position is updated as usual during the inspection, thus achieving the effect of linear velocity integration only.

[0043] It is a persistent member variable of btRigidBody, and once set, it remains unchanged throughout all four steps of the stepSimulation pipeline (force application, unconstrained prediction, constrained solution, and motion integration): in subsequent frames, regardless of the torque applied by applyTorque, the torque generated by gravitational eccentricity, or the contact impulse of the constraint solver. All modifications require factor multiplication during the motion integration phase to become effective angular velocity. A factor of zero ensures that the effective angular velocity remains constant at zero, fundamentally blocking any influence of angular velocity on attitude at the physics engine integrator level. It always remains a zero vector. This interface has been stably available since Bullet version 2.73. You can use Bullet's default btDefaultMotionState; there is no need to customize the motion state class.

[0044] Step 4: Synchronize the update layer, i.e. skip quaternion readback. After the physical step is completed, the state synchronization update is performed. For aircraft entities, the current heading angle, pitch angle and roll angle saved in the kinematic layer are directly retrieved to update the OSG scene graph, without going through the quaternion to Euler angle inverse operation.

[0045] like Figure 3 As shown, the traditional pose update path is as follows: HPR is converted to Euler angles and then to quaternions, which are input into the Bullet physics world. After Bullet integration, the quaternions are read back and then converted back to Euler angles to update the OSG scene graph. This path performs two bidirectional conversions of HPR and quaternions per frame, which poses a risk of gimbal deadlock. Loss of accuracy during readback.

[0046] The proposed approach is as follows: the kinematics layer's HPR is converted into a quaternion, input into the Bullet rigid body, and the angular velocity is frozen. Only linear velocity integration is performed. Then, the quaternion readback is skipped, and the current HPR saved by the kinematics layer is directly retrieved to update the OSG scene graph. The result is a single unidirectional conversion, no readback accuracy loss, no gimbal deadlock, and complete attitude control by the kinematics layer.

[0047] like Figure 6 As shown, in the traditional scheme The inverse kinematics formula involves matrix element extraction. The normal pose is... At that time, the yaw axis Z, pitch axis X, and roll axis Y are independent. Entering the normal decomposition branch: ; ; ; In the formula, The pitch angle, For heading angle, For roll angle, It is a rotation matrix.

[0048] When the universal joint locks up, that is... At that time, the heading axis Z coincides with the roll axis Y. It enters a deadlock branch. At this point, the heading and roll axes are coupled. Fixed to 0, only by Only the rollover value can be estimated.

[0049] At the edge of deadlock, The input value slightly exceeded Scope, via The derivative amplifies floating-point rounding errors, producing invalid values. Furthermore, when... hour, , When the denominator approaches 0, numerical noise is significantly amplified, leading to the accumulation of heading angle drift.

[0050] The advantage of this approach is that it only performs the one-way construction from HPR to quaternions, without going through... Reading back and inverse transformation fundamentally avoids the deadlock problem.

[0051] Example 1: Complete Physical Simulation and Attitude Control Process of an Aircraft Implementation scenario: An operator controls a light aircraft via a keyboard, with parameters set as follows. , , , , , , , , The W / A / D keys and up / down arrow keys control the throttle, roll, and pitch, enabling coordinated 60° turns and low-speed maneuvers. (Example provided) , , .

[0052] Step 1: Velocity separation and aerodynamic model calculation.

[0053] Read control input: Roll input -1.0 (equivalent to full A key input), pitch input 0. Current speed. , , .

[0054] Lift calculation: ,but . , , , . , .

[0055] Thrust decomposition: , .

[0056] Resistance calculation: , , .

[0057] , .

[0058] .

[0059] vertical direction balance , , Maintain level flight.

[0060] When the speed drops to 20m / s, which is below hour: , , .

[0061] The lift is only 44.4% of the weight. As it descends, the aircraft noticeably sinks.

[0062] Step 2: Kinematic layer attitude calculation.

[0063] Target roll angle ,current .

[0064] .

[0065] The temperature increments by approximately 4.8° per frame, reaching -60° in about 0.75 seconds, with no overshoot. Next frame. Updated to -4.8°, continuing to approach the target.

[0066] Automatic yaw derivation: , As roll increases, yaw rate increases.

[0067] Step 3: Kinematic layer dominant posture, i.e., unidirectional writing and freezing.

[0068] Kinematic layer calculations obtained .

[0069] Convert to radians: , , .

[0070] Construct three axial-angle quaternions: ; ; .

[0071] Quaternion multiplication .

[0072] Will Set the rotation component of the rigid body world transformation matrix, and immediately set the rigid body angular velocity to a zero vector. The aircraft completed a coordinated 60° banked turn.

[0073] Step 4: Synchronous update, i.e., skipping quaternion readback.

[0074] After the physical step, the aircraft entity directly retrieves the current heading angle of 45°, pitch angle of 5°, and roll angle of -30° stored in the kinematics layer, calls the scene graph update interface to write the position and attitude into the OSG scene graph, without going through the quaternion to Euler angle readback conversion.

[0075] If the traditional approach is used, the following will be executed: using quaternions Construct rotation matrix Extract Euler angles in the order from Z to X to Y.

[0076] .

[0077] Then, it enters the normal branch.

[0078] .

[0079] .

[0080] , a slight deviation of 0.1°.

[0081] like near , , The denominator is extremely small. and The numerical noise is significantly amplified.

[0082] This solution does not go through the above readback process; the attitude accuracy depends only on the accuracy of the double-precision quaternion construction.

[0083] Example 2: Accuracy Comparison of Traditional and Proposed Solutions in Continuous Rolling This embodiment quantitatively compares the two schemes using a complete maneuver sequence. Error accumulation effect during readback.

[0084] Maneuver sequence: The aircraft cruises at 200 m / s, initial state. The operator presses button A to perform a full 60° roll for 10 seconds, then releases the button to automatically return to the horizontal position. Then press the D key to perform a reverse full roll -60° for 10 seconds, finally returning to the starting position. Total duration: 25 seconds. Approximately 3000 frames. This example uses a small pitch angle. To avoid gimbal deadlock, focus on observing the accumulation of readback errors when the reading is well below the deadlock threshold.

[0085] Analysis of error sources in each frame readback of the traditional solution: Initial frame: Kinematic layer Construct quaternions using the Euler angle to quaternion conversion. Calculate each component according to formulas (11a) to (11c) and multiply them: , is a precise unit quaternion.

[0086] After frame 1200, i.e., a 60° roll lasting 10 seconds: the kinematic layer has... Calculate to -60° and write. This involves two rounds of conversion loops: "Euler angles to quaternions written to the physics engine, then quaternion readback, and finally Euler angle extraction." During readback, from... Extracting the rotation matrix ,calculate Matrix elements In the rotation sequence from Z to X to Y, it equals .

[0087] The theoretical value is The actual value stored in double precision is -0.0871557..., containing approximately Floating-point rounding error on the order of magnitude.

[0088] After clamping .

[0089] Enter the normal branch: ; .

[0090] This involves the problem of iterative error accumulation in double-precision floating-point operations. (Each frame) During the reread, quaternions After Bullet physical integration, each component contains approximately double-precision rounding error. Order of magnitude, quaternion rotation matrix The operation then propagates this error to each element of the matrix. numerator in calculation and denominator All contain from The propagation error, and It is sensitive to errors in both input quantities.

[0091] During the continuous rollover of 3000 frames, the tiny readback deviation of each frame accumulates frame by frame. After 1200 frames, the accumulated drift of the heading angle can reach 1 to 3°, and there is still a residual deviation of 0.5 to 1° after leveling off.

[0092] Accuracy analysis of this solution: HPR is calculated directly in each frame of the kinematic layer: Without quaternion backreading, the three variables of heading, pitch, and roll stored in the kinematics layer are directly used as update parameters for the scene graph.

[0093] After 3000 frames, the heading angle is equal to: ; In the formula The heading angle after 3000 frames. For the first The yaw increment for each frame is calculated independently using formula (10) without introducing any... Readback error. Attitude angle error originates solely from the double-precision rounding approximation of the initial HPRToQuat single-pass construction. The amplification does not accumulate during the entire 25-second maneuver. The calculated kinematic value of the heading angle after a 10-second 60° roll is completely consistent with the actual performance, and the residual deviation after leveling off is less than 0.01°, which is indistinguishable within the simulation accuracy requirements.

[0094] Those skilled in the art will understand that the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described osgEarth-based aircraft physical simulation and attitude control method.

[0095] In conjunction with the above specific implementation methods and embodiments, the osgEarth-based aircraft physics simulation and attitude control method provided by the present invention has the following significant advantages compared with the prior art: First, it completely eliminates the accuracy loss and gimbal lock risk caused by attitude back reading, and achieves extremely high-precision attitude synchronization.

[0096] As in Example 2, the traditional scheme executes in each frame. During readback, the approximation generated by integrating the quaternion through the physics engine is approximately Order-of-magnitude double-precision rounding errors will occur during the extraction of rotation matrix elements and execution. , The amplification is significant during computation, especially during sustained roll maneuvers. The function is highly sensitive to errors in both the numerator and denominator. These small deviations accumulate frame by frame, resulting in a heading angle drift of 1 to 3° after 1200 frames, with a residual deviation of 0.5 to 1° remaining after leveling off. Furthermore, when the pitch angle approaches... hour, ,lead to When the denominator approaches zero, numerical noise is amplified dramatically, even due to... The input value slightly exceeded Invalid values ​​are generated due to the range.

[0097] This invention simplifies the attitude update path from two bidirectional transformations to a single unidirectional transformation through the core mechanisms of "unidirectional writing to the kinematic layer and freezing of angular velocity" and "skipping quaternion readback". The attitude accuracy depends only on the double-precision rounding of the initial HPRToQuat single construction (approximately...). Furthermore, the error does not accumulate at all during the entire complex maneuver spanning 25 seconds and 3000 frames. Data from Example 2 conclusively demonstrates that after a 60° roll and leveling off in 10 seconds, the residual deviation of this scheme is less than 0.01°, indistinguishable within the simulation accuracy requirements. Simultaneously, because it completely eliminates... By employing inverse transformation, this invention fundamentally avoids the gimbal lock-up problem, ensuring the mathematical rigor and reliability of simulations under extreme attitudes.

[0098] Second, it significantly reduces computational complexity, perfectly meeting the demanding performance requirements of large-scale real-time situational simulation.

[0099] As shown in Example 1, this invention creatively reduces the six-degree-of-freedom rigid body dynamics of a fixed wing to horizontal velocity. and vertical velocity Two independent one-dimensional problems. By introducing stall velocity... For the segmented lift model with defined boundaries, lift calculations are simplified to pure algebraic operations. For example, in the stall region, lift is directly applied... Calculations are performed in the normal flight zone. Calculation, and supplemented by Simulates additional lift. This design completely eliminates the heavy computational burden of relying on complex aerodynamic coefficient tables and solving six-degree-of-freedom coupled differential equations, which is a traditional approach.

[0100] This dimensionality reduction and simplification process, while still accurately simulating key physical phenomena (such as in Example 1 where the speed drops to 20 m / s, the lift is accurately calculated to be 44.4% of gravity, thus correctly simulating...), allows for accurate simulation of... The system significantly reduces the computational overhead per frame by reducing the sinking acceleration. This enables the system to support real-time situational simulation of large-scale aircraft clusters with extremely low CPU utilization, solving the technical bottleneck of existing high-fidelity simulation frameworks that struggle to balance accuracy and real-time performance.

[0101] Third, it achieves smooth attitude transition without overshoot and highly stable attitude control.

[0102] This invention introduces a first-order exponential smoothing algorithm (Equation 9) and a coordinated turning mechanism (Equation 10) at the kinematic level. As shown in Example 1, at the time step... Under these conditions, the algorithm can smoothly and without overshoot approach the target roll angle of -60° in approximately 0.75 seconds, with increments of approximately 4.8° per frame. Simultaneously, through... The automatic calculation of yaw rate ensures that the flight trajectory conforms to the coordinated turning law in real aerodynamics, thus avoiding sideslip.

[0103] More importantly, by forcibly freezing the rigid body angular velocity in step three... This invention strips the physics engine of its power to generate attitude changes through torque integration, thus placing attitude control entirely under the control of the kinematic layer. This design completely isolates the interference of the physics engine's underlying numerical integration errors on the upper-level attitude control, ensuring the consistency of the aircraft's response and visual smoothness under complex and sudden control inputs, thereby enhancing the realism of the simulation.

[0104] Fourth, it improves the overall rendering consistency and engineering practicality of the digital twin system.

[0105] Traditional solutions often require forcibly fixing and estimating the heading angle when encountering gimbal lock-up edges or amplified numerical noise. This can easily lead to visual flaws in the OSG scene model, such as abrupt changes, high-frequency jitter, or clipping, severely impacting the immersive experience of the digital twin system. This invention bypasses the readback process and directly calls the scene graph update interface to write the precise HPR saved in the kinematic layer into the OSG, ensuring absolute consistency between the physical simulation data and the 3D rendered scene.

[0106] This not only eliminates visual anomalies but also significantly simplifies the data interaction interface between the two heterogeneous engines, osgEarth and Bullet. The system no longer needs to maintain a complex bidirectional state machine, reducing system integration coupling and code maintenance costs. The method provided by this invention has clear logic and a concise interface, possessing high engineering practical value and can be widely applied to digital twin systems in various fields such as civil aviation flight simulation and UAV swarm control.

[0107] 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 method for physical simulation and attitude control of aircraft based on osgEarth, characterized in that, Includes the following steps: Step 1, velocity separation and simplified aerodynamic model calculation: reduce the six-degree-of-freedom rigid body dynamics of the fixed wing to two independent one-dimensional problems of horizontal velocity and vertical velocity, and calculate the combined velocity, lift, thrust, drag and new velocity; Step 2, kinematic attitude calculation: The target attitude angle is calculated using first-order exponential smoothing, and the yaw input is derived through a coordinated turning mechanism; Step 3, One-way writing of kinematics layer and freezing of angular velocity: The target attitude angles calculated by the kinematics layer are constructed as quaternions according to the specified rotation order, written into the physics engine and the rigid body angular velocity is frozen. Step 4, Synchronous Update Layer Skips Quaternion Readback: After the physics step is completed, the current heading angle, pitch angle and roll angle saved in the kinematics layer are directly retrieved to update the OSG scene map, without going through the quaternion to Euler angle readback inverse operation.

2. The method for aircraft physics simulation and attitude control based on osgEarth according to claim 1, characterized in that, In step one, the velocity vector is synthesized. The calculation formula is: ； In the formula, The magnitude of the horizontal velocity. Vertical velocity, This is the heading angle.

3. The method for aircraft physics simulation and attitude control based on osgEarth according to claim 1, characterized in that, In step one, a segmented lift model is used, with stall speed as the starting point. As a dividing line: When horizontal speed At that time, the base lift ; when At that time, the base lift ; Additional lift ,in ; Total lift ; In the formula, For the mass of the aircraft, It is the acceleration due to gravity. The pitch angle, This is the velocity factor.

4. The method for aircraft physics simulation and attitude control based on osgEarth according to claim 1, characterized in that, In step one, the lift force is decomposed into a vertical component. and horizontal backward component ; Thrust decomposed into and ; resistance ; New speed calculation is , ; In the formula, For total lift, The pitch angle, For vertical components, The horizontal backward component, For maximum thrust, This is the throttle coefficient. For horizontal thrust, For vertical thrust, The drag coefficient, Calculate speed for resistance. As a resistance multiplier, As resistance, For a new level of speed, For horizontal acceleration, For time step, For the new vertical velocity, This is the vertical acceleration.

5. The method for aircraft physics simulation and attitude control based on osgEarth according to claim 1, characterized in that, In step two, the formula for the smooth transition of the target attitude angle is: ； The formula for coordinated turn yaw input is: ； In the formula, For the updated attitude angles, The attitude angle of the previous frame. For the target attitude angle, To approximate the rate coefficient, For time step, For yaw angle increment, The current roll angle, For the maximum roll angle, To coordinate the turning gain coefficient, This is the maximum yaw rate.

6. The method for aircraft physics simulation and attitude control based on osgEarth according to claim 1, characterized in that, In step three, the target attitude angle Quaternions are constructed by rotating from Z to X to Y. The order of quaternion multiplication is as follows: ; The components are constructed as follows: , , ; Freeze the rigid body angular velocity to ; In the formula, For the target quaternion, For the heading quaternion, For the quaternion of pitch, For rolling quaternions, For heading angle, The pitch angle, For roll angle, , , These are the radian values ​​corresponding to the heading angle, pitch angle, and roll angle, respectively. ω is the rigid body angular velocity.

7. The method for aircraft physics simulation and attitude control based on osgEarth according to claim 1, characterized in that, In step four, skipping the quaternion backread specifically means not executing the inverse formula from quaternions to Euler angles. , , ; In the formula, The pitch angle, For heading angle, For roll angle, For rotation matrix, , , , , Rotation matrices The elements in the corresponding row and column.

Images