EP4649284A1 - Verfahren und messeinrichtung zur bestimmung der lage eines objektes - Google Patents

Verfahren und messeinrichtung zur bestimmung der lage eines objektes

Info

Publication number
EP4649284A1
EP4649284A1 EP24700394.0A EP24700394A EP4649284A1 EP 4649284 A1 EP4649284 A1 EP 4649284A1 EP 24700394 A EP24700394 A EP 24700394A EP 4649284 A1 EP4649284 A1 EP 4649284A1
Authority
EP
European Patent Office
Prior art keywords
rotation
angle
cumulative
function
increments
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
EP24700394.0A
Other languages
English (en)
French (fr)
Inventor
Christoph Sieg
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Deutsches Zentrum fuer Luft und Raumfahrt eV
Original Assignee
Deutsches Zentrum fuer Luft und Raumfahrt eV
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Deutsches Zentrum fuer Luft und Raumfahrt eV filed Critical Deutsches Zentrum fuer Luft und Raumfahrt eV
Publication of EP4649284A1 publication Critical patent/EP4649284A1/de
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation

Definitions

  • Gyroscopes are used to measure the rate of rotation (angular velocity) around a predefined axis. If at least three gyroscopes are used, which are oriented so that the direction vectors along the axes of their measurement are linearly independent, the angular velocity vector ⁇ of a rotation around an axis of rotation arbitrarily oriented in space can be determined. If the gyroscopes are firmly connected to an object, such as a vehicle (e.g. a rocket), the (relative) change in position of the vehicle (orientation in space) can be calculated from the rate of rotation by integrating the signal. With JG/JG With a known initial orientation, the absolute orientation of the vehicle in space is then known.
  • a vehicle e.g. a rocket
  • the gyroscopes are usually bundled with acceleration sensors for measuring linear acceleration in a sensor package (an inertial measurement unit (IMU)), which is available as a module and is the central component for determining the orientation and position of the vehicle.
  • IMU inertial measurement unit
  • the IMU also contains signal processing stages with the aim of generating a signal at the user interface that can be used in the on-board computer. After digitizing and possibly filtering the signals, the frequency is reduced so that a sufficiently low-frequency signal that meets the specifications of the interface protocol is available for further processing at the interface.
  • this integral corresponds to the (relative) angle of rotation (angle increment) ⁇ ⁇ ⁇ , by which the IMU moves during this sampling time step ⁇ ⁇ ⁇
  • the ⁇ -th time range can be preceding, ( ⁇ -1)-th time range [ ⁇ ⁇ 2 , ⁇ ⁇ 1 ]).
  • the integral ⁇ measured by the sensors k the rotation rate ⁇ does not directly give the angle increment ⁇ ⁇ , which directly indicates the axis of rotation and the angle of rotation of the rotation that occurred in the sampling interval. Instead, the relative angle of rotation increment is ⁇ ⁇ from the integrals of the current and past rotation rates. 2.
  • the direct summation of temporally consecutive angle increments associated with the time intervals ⁇ ⁇ , ⁇ ⁇ +1 does not result which is used for rotation in the time interval [t i ⁇ 2+k ,t i+k ] angle increment. Instead of the simple summation of the angle increments requires a more complicated mathematical operation. John E. Bortz.
  • ⁇ (t) is approximated in this interval by a polynomial of (n – 1)-th degree by ⁇ ⁇ ⁇ ⁇ ⁇
  • equation (2.6) is evaluated to calculate the iterative rotation angle increment ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ This also involves increased computing effort, which increases every ⁇ -th recursion step.
  • equation (2.6) only applies to unit quaternions. Its application to non-normalized quaternions leads to an additional (numerical) error.
  • DE 102022121662 B3 discloses a method for determining the position of an object by determining the angle of rotation of the object, whereby Angle of rotation increments are measured for successive measurement points in time and combined to form a cumulative angle of rotation for a cumulative time range containing the points in time.
  • the cumulative angle of rotation for a cumulative time range is calculated as a function of exponents of an approximation of at least second order of the Baker-Campbell-Hausdorff formula for the Lie rotation group SO(3), whereby iteratively cumulative angle of rotation increments for the points in time are determined from the function of the exponents with a measured angle of rotation increment for the respective time interval and the cumulative angle of rotation increment determined for the previous point in time or a specified initial angle of rotation increment as input variables, and the cumulative angle of rotation for the time range is determined by the last cumulative angle of rotation increment determined at the last measurement point in time.
  • the iterative angle of rotation increments are calculated directly without detours via quaternions.
  • US 2018/0231385 A1 discloses a method for determining the pose parameters of an inertial measurement unit (IMU) sensor, comprising the steps of collecting measurement data generated by IMU sensors, using a processor to integrate the measurement data over time, including any errors, generating a temporally continuous error propagation model, and integrating over time.
  • the error propagation model can be used to generate compensation gradients for the pose parameters.
  • EP 1637840 A1 discloses a method for compensating conicity in a strapdown inertial navigation system that uses groups of five consecutive incremental rotation angles of a body-fixed coordinate system measured by orthogonally mounted gyros at regular measurement intervals. Each group of five measurements is obtained during a group interval corresponding to five measurement intervals.
  • the cone-compensated angular displacement of the body-fixed coordinate system around a fixed axis in space during a p-th group interval is obtained by summing the five measured incremental angles and a cone compensation term.
  • the cone compensation term consists of the sum of half the cross product of a first and a second vector sum and the weighted sum of three vector cross products.
  • the second vector sum is the sum of the five incremental angle of rotation in a group.
  • the first vector sum is the sum of the second vector sum over p groups.
  • the multiplier and the multiplicand of each vector cross product is a weighted sum of five measured incremental angles.
  • the cone-compensated angular displacement can be summed over p groups to obtain an accurate estimate of the vector angle of rotation over several group intervals.
  • the object of the present invention is to create an improved method.
  • the object is solved with the method with the features of claim 1 and the measuring device with the features of claims 11 and 12.
  • Advantageous embodiments are described in the subclaims. It is proposed that cumulative angle of rotation increments be used to determine the position of an object.
  • the rotation angle increments and rotation angles contain angle information for the required spatial directions and can have vectors with the vector components of the three units i, j, k of the unit quaternion described above or the three spatial directions x, y, z. Even if the rotation angle increments and rotation angles are not explicitly referred to as vectors in the present application and the claims for the sake of simplicity, the embodiment as a vector or representations by skew-symmetric matrices thereof is included.
  • the method has the advantage that the iterative rotation angle increments are calculated directly without detours via quaternions.
  • absolute angle increments and thus in particular the angle of rotation parameterizing the position of the object (body) can also be calculated with a given precision.
  • the only transcendental function to be calculated is ⁇ ( ⁇ ) is nevertheless approximated with high precision by a rational function and can therefore be calculated efficiently and with a previously known number of calculations.
  • the translational speed and position can be calculated from the accelerations/relative speed increments measured by the IMU's acceleration sensors. In combination with the position and a model for gravity, complete inertial navigation within the IMU appears to be possible, independent of the on-board computer.
  • Trigonometric functions and rotation matrices can also be calculated efficiently and precisely in other application areas such as computer graphics/computer vision, CAD for given angles of rotation.
  • the calculation of the cumulative angle of rotation increments ⁇ ⁇ for a particular measurement point ⁇ ⁇ can be considered as a function of at least second order in the second Argument approximated exponents of the Baker-Campbell-Hausdorff formula for the Lie rotation group SO(3).
  • Calculating the function ⁇ in line 4 of Algorithm 1 requires in particular the approximate calculation of only a single transcendental function ( 1.3) which is the value of the last cumulative angle increment, i.e.
  • the algorithm is based on a simpler iteration and should simplify the analysis of the propagation of measurement uncertainties and measurement noise compared to the quaternion-based algorithm described in the introduction to the state of the art. These properties allow a direct implementation on the hardware (for example a Field Programmable Gate Array (FPGA)) of the IMU.
  • the algorithm can thus be implemented as part of a sensor package and the signal processing of the sensor data for output to the user interface.
  • the cumulative angle of rotation ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ for the time range ⁇ is a function of the last cumulative angle increment ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ , which for the last time ⁇ ⁇ ⁇ of the interval of the time range I. It can be assigned to this cumulative angle increment ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ ⁇
  • the integral ⁇ measured by the sensors k the rotation rate ⁇ does not directly give the angle increment ⁇ the direct axis of rotation and the angle of rotation of the ⁇ rotation that occurred in the sampling interval. Instead, ⁇ from the ⁇ To calculate integrals of the current and past rotation rates. 2.
  • the direct summation of temporally consecutive time intervals ⁇ ] associated angle increments ⁇ ⁇ , ⁇ ⁇ does not yield the rotation in the time interval [ ⁇ meas, ⁇ 2, ⁇ meas, ⁇ ] belonging to the angle increment.
  • a more complicated mathematical operation must be carried out.
  • the method is based on the complex mathematical operations required for time-varying axes of rotation instead of a simple summation of the angle increments.
  • the reason for this is that rotations around non-parallel axes of rotation do not interchange with each other. Rather, the rotations in (three-dimensional) space form a group, the rotation group SO(3), which is a so-called Lie group.
  • the mathematical area of the Lie groups forms the theoretical framework from which the representations of the group elements in the form of rotation matrices or in the form of unit quaternions arise.
  • the rotation matrices are special (their determinant is 1) orthogonal (their transpose is their own inverse) matrices in three dimensions.
  • the name of the group, SO(3) is derived from the terms special orthogonal rotation matrices in 3 dimensions.
  • Unit quaternions are the extension of complex numbers to three independent imaginary units, whose norm calculated from the real and imaginary parts is 1 and are multiplied by means of a non-commutative, norm-preserving product.
  • a rotation can be represented by a unit quaternion, whereby the representation is not unique, but both q and -q can be used to represent a specific rotation.
  • the operations described above for calculating the angle increments and their summation then follow.
  • Number n of rotation angle increments ⁇ ⁇ , which were measured at consecutive times tmeas,k with k 1, ..., n, and the cumulative angle increment determined for the previous time interval ⁇ ⁇ 1
  • the calculation is based on only one transcendental function, which reduces the time and computational effort.
  • Figure 2 for the position of an object
  • Figure 4 Diagrams of the relative errors when using the Padé approximation to calculate the transcendental function ⁇ ( ⁇ )
  • Figure 5 Diagrams of the relative errors when using the Padé approximation directly from ⁇ ( ⁇ ) and a regularized expression to calculate the transcendental function ⁇
  • Figure 6 Relative error of approximation of ⁇ ⁇ ⁇ and ⁇ ⁇ ⁇ ⁇ 0 ( ⁇ ) which on a Padé approximation of ⁇ ( ⁇ )
  • Figure 7 Relative error of the
  • FIG. 2 shows a sketch of a measuring device 1 for determining the position of an object 2.
  • the measuring device has at least one inertial measuring unit IMU, which is used to measure the angle of rotation increments ⁇ ⁇ ⁇ consecutive points in time ⁇ ⁇ for the detection of rotating (circular) movements of the inertial measuring unit IMU and the associated object 2 in the three mutually orthogonal spatial axes (X or Y or Z axis) of the Cartesian coordinate system is set up.
  • ⁇ It is optionally conceivable to measure the accelerations or relative speed increments as components along the three spatial axes for the speed and, if applicable, position and to feed them to the data processing unit.
  • the data processing unit 3 can, for example, be a microprocessor or microcontroller that executes a computer program with program commands that cause the microprocessor or microcontroller to execute the method steps of the method according to the claim.
  • the data processing unit 3 can, however, also be designed as hard-wired hardware logic, such as a field programmable gate array or similar.
  • the method enables a direct calculation of the cumulative angle of rotation ⁇ ⁇ iteratively from the angle of rotation increment saved in the last step and the currently measured one ⁇ ⁇ ⁇ carried out.
  • the rotation matrix obtained in this process can be used to calculate the speed and, if necessary, position from the accelerations or speed increments that are also measured.
  • BCH Baker-Campbell-Hausdorff
  • ⁇ ( ⁇ , ⁇ ) is given by the Baker-Campbell-Hausdorff formula as a formal series expansion in the total order m of ⁇ and ⁇ and therefore has the following schematic form: ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ with the coefficient ⁇ ( ⁇ ⁇ , ⁇ ⁇ ) , where the actual form is due to the general non-commutability of the group elements ⁇ ⁇ ⁇ ⁇ ⁇ is much more complicated.
  • the general Baker-Campbell-Hausdorff (BCH) formula is not yet available as a formal series in a form that is suitable for the summation of temporally consecutive series with the time intervals [t i ⁇ 3+k , t i-2+k ], [t i-2+k , t i-1+k ] associated angle increments ⁇ ⁇ 1 , ⁇ ⁇ ⁇ is suitable to determine the actual rotation in the time interval [ti ⁇ 3+k, ti-1+k] belonging to the rotation angle increment, i.e. the rotation angle ⁇ ⁇ ⁇ O ⁇ ut for the time range I (in this case consisting of two subintervals).
  • the BCH formula still has to be specialized for the special case of the Lie group SO(3).
  • a suitable evaluation / calculation of the series expansion must be carried out.
  • SO(3) a suitable evaluation of the expression, which is still present as a series, is carried out.
  • the specialization to the Lie group SO(3) initially only provides the required relation for the summation of two angle elements, for the previously described case ⁇ ⁇ where ⁇ is still given as a series expansion in the total order m as in equation (3.2) and ⁇ ⁇ is the newly cumulated rotation angle increment resulting from the newly measured relative rotation angle increment ⁇ ⁇ and the previously calculated cumulative rotation angle increment.
  • the only function contained in the system of equations and to be calculated by approximation is the transcendental function ⁇ ( ⁇ ) .
  • the expressions are whereas in the first argument ⁇ ⁇ exactly.
  • Figure 2 shows a corresponding flow chart of this method or the above procedure for determining a cumulative rotation angle increment from rotation angle increments measured with a rotation angle sensor for successive time intervals.
  • the starting point is the exponential mapping, which is also defined for arguments ⁇ , the elements of a Lie algebra, i.e. the tangent space of a Lie group at Element of the identity and which can be represented by matrices.
  • the exponential mapping is defined by its Taylor series ⁇ ⁇ where ⁇ is an element of a representation of a Lie algebra. This element ⁇ can be given, for example, by the pure quaternions in the case of rotations. Powers in the above equation are to be understood as being carried out with the corresponding product, which is defined for the representation, e.g. the matrix product in the case of a matrix representation or the quaternion product in the case of a representation by quaternions.
  • ⁇ ⁇ ( ⁇ , ⁇ ) ⁇ ⁇ ⁇ ⁇ .
  • the exponent ⁇ is expanded in orders of ⁇ , as is shown, for example, in M. Müger's "Note on the theorem of Baker-Campbell-Hausdorff-Dynkin".
  • Explicit expression for the BCH formula up to the second order, m ⁇ 0, 1, 2
  • the right fist or corkscrew rule can be used here, where the thumb of the right hand points in the direction of the vector and the fingertips of the clenched fist indicate the direction of rotation for rotations with a positive angle, i.e. the norm of the vector-valued angle of rotation.
  • ⁇ ′ ⁇ ′ ⁇ ⁇ ′ parameterized, whose norm 0 ⁇ ⁇ ′ ⁇ ⁇
  • ⁇ ⁇ ′ ( ⁇ ⁇ 2 ⁇ 0) ⁇ ⁇ (3.1)
  • , ⁇ ⁇ ′ sgn( ⁇ ⁇ 2 ⁇ 0 ) ⁇ ⁇ (3.2) given, i.e.
  • the representative of minimal norm can then be constructed by transforming equation (3.5) as follows —
  • integer values ⁇ ’ over the interval 1 ⁇ ⁇ ′ can be iterated, either starting forward from the bottom or ⁇ from the upper interval limit.
  • the current value can be ⁇ ' respectively calculated and the sign change depends on the size ⁇ ′
  • the largest value of ⁇ ', for the ⁇ ′ ⁇ 0 is valid, then the value for ⁇ and ⁇ ⁇ ⁇ 0.
  • the representative with the minimum norm can then be calculated using equation (3.10).
  • the expression (3.10) for calculating the representative with the minimum norm requires the calculation of the square root ⁇ 1 + ⁇ .
  • ⁇ ( ⁇ ) is at any time t k , to which a new relative angle increment ⁇ is available for ⁇ the calculation of the cumulative angle increment ⁇ ⁇ required.
  • the application of the algorithm described above for renormalization of the angle increments guarantees that the norm ⁇ ⁇ 1 the condition 0 ⁇ ⁇ ⁇ ⁇
  • numerical approximations for ⁇ ( ⁇ ) which allow ⁇ ( ⁇ ) efficiently and with the desired precision over the entire required range of values.
  • Equation (1.3) has the following Taylor series expansion: evaluated ⁇ N express.
  • the Riemannian is defined over the series (3.21) Inserting equation (3.20) into equation (3.19) then gives for the Taylor series from ⁇ ( ⁇ ) where ⁇ ⁇ 0+1 the remainder of the approximation of ⁇ by the Taylor polynomial of degree 2 ⁇ 0 in ⁇ It is important that ⁇ (2n), i.e. the Riemann not with the argument ⁇ is confused.
  • the Padé approximation proposed by the invention is considerably more accurate than using the Taylor polynomial for the approximation, with approximately the same amount of computational effort.
  • the function and the numerator polynomial and the denominator polynomial have the degrees ⁇ 0 +1 or1 in ⁇ 2
  • Figure 5 shows the relative errors when directly using the Padé approximation (3.25) to calculate ⁇ (2) ( ⁇ ) and when using the regularized approximation ⁇ ( ⁇ 2 ⁇ ) ⁇ , ⁇ 0( ⁇ ) from equation (3.40) for different values of ⁇ 0 over the entire relevant value range 0 ⁇ ⁇ shown.
  • Diagram a) shows the relative errors where the Using the approximation of ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ 0 from equation (3.34) into the Equation (3.36) is given.
  • Diagram b) shows the relative error when given by the regularized approximation from equation (3.40) is.
  • the accuracy is determined by the Termination of the expansion of equation (3.41) according to the second order in ⁇ 2 ⁇ ⁇ ⁇ limited, which is the most relevant case in practice ⁇ ⁇ ⁇ 1 However, it does cover this. If necessary, an extension of equation (3.41) can be made using the methods described above without any conceptual obstacles.
  • the calculation of trigonometric functions and the rotation matrix is explained below.
  • the acceleration sensors in the IMU sensor package either measure the acceleration directly or they measure relative speed increments, which are given in analogy to relation (1.1) for the relative angle increments by are given.
  • the acceleration of the acceleration sensors with respect to an inertial system I expressed in the reference system S of the sensor at the same time t at which the value of the acceleration a(t) is present.
  • equation (1.1) therefore also includes rotational movements that the sensor performs during the integration time.
  • the speed increments according to equation (1.1) cannot therefore be used directly to determine the speed and position of the vehicle with respect to an inertial system or other reference system. Instead, speed increments are required that have been calculated by integrating the acceleration in a coordinate system that is fixed at least in the integration interval.
  • the individual speed increments are again integrated into a common reference system, e.g. S I-1 at the start time ⁇ the summation, to transform.
  • the summation is (without gravity)
  • ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ 1 1+ ⁇ ⁇ which is the starting time ⁇ ⁇ 1 from cumulative relative Angle increment ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ 0 ⁇ ⁇ associated rotation matrix.
  • the algorithm is not only relevant for the application described above, the processing of the sensor data of an IMU and the efficient calculation of the position and attitude of a vehicle.
  • the efficient calculation of trigonometric functions presented here should have a wide range of applications in various areas, e.g. in computer vision. It is shown below that the algorithm used for updating the cumulative angle increment ⁇ ⁇ already calculated result for the function value of the function ⁇ ( ⁇ ) can be used to efficiently and especially without approximating any other transcendental functions the (trigonometric) functions ⁇ ⁇ 2 ⁇ ⁇ 2 ⁇ and ⁇ ⁇ ⁇ 2 ⁇ .
  • the associated rotation matrix ⁇ ⁇ ⁇ ⁇ ( ⁇ ) of the active rotation between the coordinate systems F and G are mapped onto each other by means of rotation.
  • Knowledge of this matrix is an essential building block for calculating the translational velocity and ultimately the position from the accelerations or relative velocity increments measured by the acceleration sensors.
  • Figure 6 shows the relative errors of the approximations of ⁇ (Diagram a)) and cos ⁇ (Diagram b)) with the Padé approximation according to equation (3.34) in the approximations according to equation (3.46) in the relevant value range 0 ⁇ ⁇ ⁇ ⁇ .
  • Figure 7 the relative errors of the approximations of sin ⁇ (diagram a)) and cos ⁇ (diagram b)) in the direct approximation by Taylor polynomials of comparable degree.
  • the use of the Padé approximation by equation (3.34) in the approximations derived in this section according to equation (3.46) provides a more accurate approximation over the entire range of values than a direct approximation by a Taylor polynomial of comparable degree.
  • Algorithm 6 Calculation of ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ out of ⁇ , ⁇ 2 and B ⁇ 2 ⁇ ⁇ 1: procedure ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ 2 , B ⁇ 2 ⁇ 2: Calculate ⁇ 3: Calculate ⁇ ⁇ 2 ⁇ out of 2 ⁇ ⁇ ⁇ 4: Calculate Calculate the rotation matrix ⁇ 2 ⁇ ⁇ ⁇ , 2 2 5: ⁇ 2 ⁇ ⁇ ⁇ , ⁇ ⁇ 2 ⁇ using equation (3.50)
  • Algorithm 7 is a slight modification of Algorithm 6, which shows how the coordinates ⁇ ⁇ of a vector ⁇ based on equation (3.51) directly in the reference system G from the coordinates ⁇ ⁇ in the
  • Algorithm 7 Calculation of ⁇ ⁇ out of ⁇ , ⁇ 2 and B ⁇ 2 ⁇ , ⁇ ⁇ 4: Calculate calculate ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ 2 5: Be ⁇ ⁇ ⁇ , ⁇ , ⁇ 2 ⁇ , ⁇ ⁇ using equation (3.51)
  • Both 7, e.g. can be used for the efficient and precise calculation of trigonometric functions and rotations in real-time applications, e.g. for calculating rotations in the field of computer graphics / computer vision and CAD.
  • Figure 9 shows a flow chart of the method in which at least one gyroscope measures the angle increments in the x-direction, in the y-direction and in the z-direction.
  • the three measured values for the angle increments ⁇ ⁇ 1 ⁇ ,2,3 can be vectorized to create a rotation angle increment vector ⁇ ⁇ This allows the individual angle increments ⁇ ⁇ recursively from the exponent ⁇ ⁇ ⁇ ⁇ ⁇ with the delay z -1 for inclusion of the previous rotation angle increment vector ⁇ can be calculated.
  • the last cumulative angle of rotation increment calculated in this way then forms the cumulative angle of rotation ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ as a measurement result.

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Complex Calculations (AREA)

Abstract

Es wird ein Verfahren zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehwinkel des Objektes beschrieben, wobei relative Drehwinkelinkremente Δθ k an aufeinanderfolgenden Messzeitpunkten t k mit k = 1,..., n gemessen und für die jeweiligen Messzeitpunkte t k zu kumulierten Drehwinkelinkrementen θ k zusammen- gefasst werden. Es erfolgt ein Berechnen der kumulierten Drehwinkelinkremente θ k für die Messzeitpunkte t k iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen Δθ k , die an aufeinanderfolgenden Messzeitpunkten t k mit k = 1,..., n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θ k-1 und der transzendenten Funktion (I) mit ϛ als Funktion des kumulierten Drehwinkelinkrementes θ k-1 zum jeweils vorhergehenden Messzeitpunkt t k -1, wobei ein Approximieren der transzendenten Funktion (I) durch Padé-Approximationen. Durch diese Padé-Approximationen der transzendenten Funktion kann auch sehr ressourceneffizient und genau eine automatische rechnergestützte Berechnung von trigonometrischen Funktionen von Winkeln θ k oder Drehungen aus gegebenen Winkeln θ k vorgenommen werden.

Description

Deutsches Zentrum für Luft- und Raumfahrt e.V. Anwaltsakte: Königswinterer Straße 522-524 V/DLR-0799-WO 53227 Bonn-Oberkassel 2022/224 Deutschland Datum: 8. Januar 2024 Verfahren und Messeinrichtung zur Bestimmung der Lage eines Objektes Das Projekt, das zu dieser Anmeldung geführt hat, wurde durch das Forschungs- und Innovationsprogramm Horizont 2020 der Europäischen Union unter der Fördervereinbarung Nr.101004205 gefördert. Die Erfindung betrifft ein Verfahren zum Bestimmen der Lage eines Objektes aus den Drehraten des Objektes, wobei relative Drehwinkelinkremente ^^ an ^^ aufeinanderfolgenden Messzeitpunkten ^^ ^^ mit k = 1, …, n gemessen und für die jeweiligen Messzeitpunkte ^^ ^ zu kumulierten Drehwinkelinkrementen ^^ zusamm ^ ^^ en- gefasst werden. Die Erfindung betrifft weiterhin eine Messeinrichtung zum Bestimmen der Lage eines Objektes aus den Drehraten des Objektes, wobei die Messeinrichtung eine Sensoreinheit zum Messen von relativen Drehwinkelinkrementen ^^ ^^ für aufein- anderfolgende Messzeitpunkte ^^ ^ mit k = 1, …, n und eine Datenverarbeitungseinheit ^ hat, die eingerichtet ist, die an aufeinanderfolgenden Messzeitpunkten ^^ ^^ gemessenen relativen Drehwinkelinkremente ^^ ^^ für die jeweiligen Messzeitpunkte ^^ ^^ zu kumulierten Drehwinkelinkrementen ^^ ^^ zusammenzufassen. Gyroskope dienen der Messung der Drehrate (Winkelgeschwindigkeit) um eine vordefinierte Achse. Bei Verwendung mindestens dreier Gyroskope, die so orientiert sind, dass die Richtungsvektoren entlang der Achsen ihrer Messung linear unabhängig sind, lässt sich der Winkelgeschwindigkeitsvektor ω einer Drehung um eine beliebig im Raum orientierte Drehachse bestimmen. Wenn die Gyroskope fest mit einem Objekt, wie insb. einem Fahrzeug (z. B. einer Rakete) verbunden sind, lässt sich aus der Drehrate durch Integration des Signals die (relative) Lageänderung des Fahrzeugs (Orientierung im Raum) berechnen. Bei JG/JG bekannter Anfangsorientierung ist dann die absolute Orientierung des Fahrzeugs im Raum bekannt. Die Kenntnis der Orientierung ist die Eingangsgröße für die Lage- regelung und für die Bestimmung der Position, wenn orientierungsabhängige Kräfte auf das Fahrzeug wirken, z. B. verursacht durch fahrzeugfeste Antriebs- oder Bremssysteme. Üblicherweise werden die Gyroskope mit Beschleunigungssensoren zur Messung der Linearbeschleunigung in einem Sensorpaket (einer Inertialmesseinheit bzw. englisch: Inertial Measurement Unit (IMU)) gebündelt, welches als Modul verfügbar ist und zentraler Baustein für die Bestimmung der Orientierung und Position des Fahrzeugs ist. Die IMU enthält neben den Sensoren auch noch Stufen der Signal- verarbeitung, mit dem Ziel an der Benutzerschnittstelle ein im Boardcomputer weiterzuverwendendes Signal zu erzeugen. Nach der Digitalisierung und evtl. Filterung der Signale erfolgt eine Reduktion der Frequenz, so dass an der Schnittstelle ein die Spezifikationen des Schnittstellenprotokolls erfüllendes hinreichend niederfrequentes Signal zur Weiterverarbeitung zur Verfügung steht. Im Fall der Gyroskope werden oft nicht direkt die Drehraten ω gemessen, sondern ihre Integrale, d. h. die über einen Abtastzeitschritt Δ ^^ ^^ integrierte (d. h. kumulierte) Drehrate, die z. B. im Zeitschritt von ti−1 nach ti, i ∈ N durch ^ ^^ ^^ gegeben ist. Im Fall einer zeitlich konstanten Drehachse entspricht dieses Integral dem (relativen) Drehwinkel (Winkelinkrement) ^^ ^^ ^^, um den sich die IMU während dieses Abtastzeitschrittes Δ ^^ ^^ gedreht hat. Soll nun an der Benutzerschnittstelle ein kumulierter Drehwinkel Δ ^^ ^ ^ ^ ^ ^^ ^^ ausgegeben werden, der sich auf einen längeren Zeitbereich ∆ ^^ ^^, ^^ ^^ ^^ = ∑ ^ ^^ ^ =1 Δ ^^ ^^ , n ∈ N von nach ^ (also eine gegenüber der Eingangsfrequenz ^^ niedrigere Frequenz < bezieht, so sind die n individuellen Drehwinkel ∆ ^^ ^^―1+ ^^ , k = 1,..., n, die innerhalb des Zeitintervalls [tI−1, tI ] gemessen werden, aufzuaddieren. Zum Verständnis der Notation wird darauf hingewiesen, dass das Zeitintervall [ ^^ ^^―1, ^^ ^^], das im Folgenden auch als Zeitbereich I bezeichnet wird, durch die Folge von ^^ Messzeitpunkten ^^meas, ^^, ^^ = ^^ wobei ^^ ^^ = ^^ ^^ ist, in ^^ Zeitintervalle [ ^^meas, ^^―1, ^^meas, ^^] aufgeteilt wird, wobei die individuellen Längen Δ ^^ ^^ = ^^meas, ^^ ― ^^ unterschiedlich sein können. Der ^^ -te Zeitbereich kann sich an vorausgehenden, ( ^^-1)-ten Zeitbereich [ ^^ ^^―2, ^^ ^^―1]) anschließen. Beträgt beispielsweise die Abtastfrequenz der Gyroskope f = 2 kHz, so dass die gemessenen Winkelinkremente auf eine konstante Abtastrate von Δt = f -1 = 5 * 10−4 s bezogen sind und soll die Frequenz an der Benutzerschnittstelle fout = 100 Hz betragen, so sind jeweils ^^ = ^^ ^^ ^^ ^^ ^^ = 20 der mit Zeitabständen Δt = 5 * 10−4 s gemessenen Drehwinkelinkremente ∆ ^^ ^^ zu summieren, um das Winkelinkrement bezogen auf einen Zeitschritt Δtout = 1 / fout = 0,01 s am Ausgang zu erhalten, also den relativen Drehwinkel, um den sich die IMU während eines Zeitschritts Δtout gedreht hat. Diese Überlegungen sind allerdings nur korrekt, wenn die Rotation um eine zeitlich konstante Drehachse erfolgt. Diese Annahme ist jedoch in der Regel nicht erfüllt. Rotiert beispielsweise ein starrer Körper ohne Einwirkung äußerer Drehmomente um eine Achse, die keine der drei Hauptträgheitsachsen ist, so führt der Winkel- geschwindigkeitsvektor im körperfesten Bezugssystem eine Präzessionsbewegung aus und seine Richtung ist nicht zeitlich konstant. Im Kontext der Navigation wird dieses Phänomen als „Coning“ bezeichnet. Im Fall einer zeitlich veränderlichen Drehachse ergeben sich nun folgende Änderungen: 1. Das von den Sensoren gemessene Integral Δφk der Drehrate ω ergibt nicht direkt das Winkelinkrement ∆ ^^ ^^, das direkt Drehachse und den Drehwinkel der im Abtastintervall erfolgten Drehung angibt. Stattdessen ist das relative Drehwinkelinkrement ∆ ^^ ^^ aus den Integralen der aktuellen und vergangenen Drehraten zu berechnen. 2. Die direkte Summation zeitlich aufeinanderfolgender mit den Zeitintervallen assoziierter Winkelinkremente ∆ ^^ ^^ , ∆ ^^+1 ergibt nicht das zur Drehung im Zeitintervall [ti−2+k,ti+k] gehörende Winkelinkrement. Anstelle der einfachen Summation der Winkelinkremente ist eine kompliziertere mathematische Operation auszuführen. John E. Bortz. „A New Mathematical Formulation for Strapdown Inertial Navigation“, In: IEEE Transactions on Aerospace and Electronic Systems AES-7.1 (1971), pp. 61–66. DOI: 10.1109/TAES.1971.310252, schlägt eine (nichtlineare) Differential- gleichung für den Drehwinkel θ(t) als Funktion der Zeit vor, die von der Drehrate ω(t) abhängt und wie folgt lautet: (2.1) In der obigen Gleichung ist der mehrkomponentige Drehwinkel θ(t) in Fettdruck dargestellt, um ihn von seiner Norm θ = ‖ ^^‖2 zu unterscheiden. Diese Differentialgleichung kann für kleine Drehwinkel θ vereinfacht werden, indem sie um Drehwinkel θ = 0 in eine Taylorreihe entwickelt wird, wobei nur die führende(n) Ordnung(en) beibehalten werden. Sie lässt sich dann beispielsweise in einem Zeitintervall [tI-1, tI], also für das Winkelinkrement ^^ ^^ ^ ^ ^ ^ ^^ ^^ lösen. Hierzu wird ω(t) in diesem Intervall durch ein Polynom (n – 1)-ten Grades approximiert durch ^^ ^^ ^ Dessen n Koeffizienten können aus n gemessenen integrierten Drehraten Δφk, k = 1, ... , n bestimmt werden. Hierzu müssen die Koeffizienten in die Polynom- approximation eingesetzt, das Integral berechnet und dann die ^^ = 1,..., ^^ Gleichungen nach den Koeffizienten aufgelöst werden. Dieses Verfahren benötigt bereits n gemessene Winkelinkremente Δφk von Zeitintervallen, die zusammen den Zeitbereich des kumulierten Drehwinkels ∆ ^^ ^^ bilden, und liefert den kumulierten Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ des Intervalls [tI−1, tI]. M. B. Ignagni. „Efficient class of optimized coning compensation algorithms“, In: Journal of Guidance, Control, and Dynamics 19.2 (1996), pp.424–429. DOI: 10.2514/3.21635. eprint: https://doi.org/10.2514/3.21635 beschreibt einen Algorithmus höherer Ordnung. Allerdings steigt für höhere Ordnungen ^^ der Rechenaufwand stark an und die Approximation der Messdaten durch Polynome hohen Grades ^^ ― 1 führt zu starken Oszillationen an den Rändern des Interpolationsintervalls, die als Runges Phänomen bekannt sind. R. A. McKern. „A study of transformation algorithms for use in a digital computer“, MA thesis. Massachusetts Institute of Technology, Dept. of Aeronautics and Astronautics, 1968. URL: http://hdl.handle.net/1721.1/14164 offenbart einen Algorithmus zur Kombination der relativen Rotationen aus aufeinanderfolgenden Zeitintervallen. Dieser Algorithmus berechnet über den Umweg der Einheits- quaternionen aus den Winkelinkrementen zweier zeitlich aufeinanderfolgender Zeitintervalle das Einheitsquaternion der relativen Lageänderung während des gesamten Zeitintervall. Dabei wird der zuvor beschriebene Algorithmus zur Berechnung des kumulierten Drehwinkels ^^ ^^ ^ für das Intervall tI] genutzt, allerdings nur für n = 2 aufeinanderfolgende Winkelinkemente. Dann wird mit dem erhaltenen Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ das Einheitsquaternion der relativen Lageänderung in dem Zeitintervall tI] berechnet zu: ^ Der Drehwinkel-Vektor im Argument der Exponentialfunktion und auf der rechten Seite der Gleichung ist als sogenanntes reines Quaternion zu interpretieren, das nur einen mit den Vektorkomponenten des Winkelinkrements ∆ ^^ ^ ^ ^ ^ ^^ ^^ gebildeten Imaginärteil besitzt. Darüber hinaus werden die trigonometrischen Funktionen durch ihre Taylorpolynome bis zur 3. Ordnung approximiert. Das Ergebnis wird dann verwendet, um ausgehend von einem Startzeitpunkt ^^0 und dem Quaternion 0 0 ^^ = 1 das Quaternion der relativen Lageänderung rekursiv zu berechnen gemäß: ^^ ^^ ^ Hierbei sind die Quaternionen mittels des Quaternionproduktes zu multiplizieren. In Analogie zu der Formel (2.3) werden die zu den Winkelinkrementen ^^ ^^ ^^―1,  ^^ ^^ ^^ gehörenden Quaternionen ^ ^^ ^^ berechnet und dann in Analogie zur Formel ^^―1 (2.4) multipliziert, um das ^ ^^^―2 ^^ = ^ ^ ^ ^ ― ― 1 2 ^^ ^ ^ ^ ^―1 ^^. (2.5) zu erhalten. Anschließend ist dann das Winkelinkrement der kombinierten Drehungen durch Inversion der Funktion (2.3) zurückzugewinnen. Die entsprechende Relation lautet für das aus den Vektorkomponenten des Winkelinkrements ^^ ^^ ^ ^ ^ ^ ^^ ^^ gebildete Quaternion wie folgt: Diese Lösung für den Spezialfall von n=2 zeitlich aufeinanderfolgenden Winkel- inkrementen ^^ ^^ ^^―1,  ^^ ^^ ^^kann auf beliebige Anzahl n von Winkelinkrementen für einen Zeitbereich und einen zugehörigen Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ verallgemeinert werden. Die iterative Berechnung des Drehwinkels ^^ ^^ ^ ^^^ = ^^ des ^^ ^^ ^ ^^ ^^ ^^ ^^ ^^ durch die Prozedur schrieben Setze ^^ Für ^^ = Berechne ^ ^^ aus ^^ ^^ ^^ mit der Gleichung (2.3) Be 0 ^^ ^^ und ^^ mit der Gleichung (2.5) Berechne ^^ aus ^ 0 ^ ^^ mit der Gleichung (2.6) Die Berechnung des Drehwinkelinkrementes ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ erfolgt indirekt über Quaternionen und beruht daher auf den drei Gleichungen (2.3), (2.4) und (2.5). Dies erfordert die Berechnung verschiedener transzendenter Funktionen (cos, sin, arcsin) z. B. nach Approximation durch ihre Taylorpolynome, wodurch der Rechenaufwand erhöht wird. Nach jeder ^^-ten rekursiven Aktualisierung des Quaternions wird die Gleichung (2.6) ausgewertet, um aus dem Quaternion das iterative Drehwinkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ zurückzugewinnen. Auch dies ist mit erhöhtem Rechenaufwand verbunden, der jeden ^^-ten Rekursionsschritt anfällt. Bei Approximation der trigonometrischen Funktionen in der Gleichung (2.3) ist das Quaternion kein Einheitsquaternion mehr. Die Relation der Gleichung (2.6) gilt jedoch nur für Einheitsquaternionen. Ihre Anwendung auf nicht normierte Quaternionen führt zu einem zusätzlichen (numerischen) Fehler. Die Berechnung der Gleichung (2.6) bei hinreichend großen Winkelinkrementen ^^ ^^ ^ ^ ^ ^ ^^ ^^, wie sie z. B. entstehen, wenn eine große Anzahl n von Inkrementen ^^ ^^ ^^ summiert wird, erfordert bei Approximation durch ihr Taylorpolynom eine Berücksichtigung recht hoher Ordnungen, um die gewünschte numerische Genauigkeit zu erreichen. Dieses erhöht den Rechenaufwand. Die bekannten Verfahren zur Bestimmung der Lage aus einem Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ eines Zeitbereichs, bspw. des Intervalls , tI] durch Zusammenfassung einer Anzahl von n Winkelinkrementen Δφ k an Zeitpunkten ^^ ^^, bspw. ^^ = die zusammen den Zeitbereich bilden, sind sehr rechenaufwändig und benötigen daher relativ viel Rechenzeit, Rechenkapazität und erfordern eine relativ leistungsfähige und damit aufwändige Elektronik mit Hardware-Rechenleistung. DE 102022121662 B3 offenbart ein Verfahren zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehwinkel des Objektes, wobei Drehwinkelinkremente für aufeinanderfolgende Messzeitpunkte gemessen und zu einem kumulierten Drehwinkel für einen die Zeitpunkte beinhaltenden Kumulations- Zeitbereich zusammengefasst werden. Es erfolgt ein Berechnen des kumulierten Drehwinkels für einen Kumulations-Zeitbereich als Funktion von Exponenten einer Approximation mindestens zweiter Ordnung der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3), wobei iterativ kumulierte Drehwinkelinkremente für die Zeitpunkte jeweils aus der Funktion der Exponenten mit jeweils einem gemessenen Drehwinkelinkrement für das jeweilige Zeitintervall und dem für den vorhergehenden Zeitpunkt bestimmten kumulierten Drehwinkelinkrement oder einem vorgegebenen Ausgangs- Drehwinkelinkrement als Eingangsgrößen bestimmt werden und der kumulierte Drehwinkel für den Zeitbereich durch das am letzten Messzeitpunkt bestimmte letzte kumulierte Drehwinkelinkrement bestimmt ist. Die Berechnung der iterativen Drehwinkelinkremente erfolgt direkt ohne Umwege über Quaternionen. US 2018/0231385 A1 offenbart ein Verfahren zum Bestimmen der Posenparameter eines Inertialmesseinheitssensors (IMU) mit den Schritten des Sammelns von durch IMU-Sensoren erzeugten Messdaten, der Verwendung eines Prozessors zum zeitlichen Integrieren der Messdaten, einschließlich etwaiger Fehler, des Erzeugens eines zeitlich kontinuierlichen Fehlerausbreitungsmodells und des zeitlichen Integrierens. Mit dem Fehlerausbreitungsmodell lassen sich Kompensations- gradienten für die Posenparameter generieren. EP 1637840 A1 offenbart ein Verfahren zur Kompensation der Konizität in einem Strapdown-Inertialnavigationssystem, das Gruppen von fünf aufeinanderfolgenden inkrementellen Drehwinkeln eines körperfesten Koordinatensystems verwendet, die von orthogonal montierten Kreiseln in regelmäßigen Messintervallen gemessen werden. Jede Gruppe von fünf Messungen wird während eines Gruppenintervalls erhalten, das fünf Messintervallen entspricht. Die kegelkompensierte Winkelver- schiebung des körperfesten Koordinatensystems um eine feste Achse im Raum während eines p-ten Gruppenintervalls wird durch Summieren der fünf gemessenen Inkrementalwinkel und eines Kegelkompensationsterms erhalten. Der Konus- kompensationsterm besteht aus der Summe von der Hälfte des Kreuzprodukts einer ersten und einer zweiten Vektorsumme und der gewichteten Summe von drei Vektorkreuzprodukten. Die zweite Vektorsumme ist die Summe der fünf inkrementellen Drehwinkel in einer Gruppe. Die erste Vektorsumme ist die Summe der zweiten Vektorsumme über p-Gruppen. Der Multiplikator und der Multiplikand jedes Vektorkreuzprodukts ist eine gewichtete Summe von fünf gemessenen Inkrementalwinkeln. Die kegelkompensierte Winkelverschiebung kann über p- Gruppen summiert werden, um eine genaue Schätzung des Vektordrehwinkels über mehrere Gruppenintervalle zu erhalten. Ausgehend hiervon ist es Aufgabe der vorliegenden Erfindung, ein verbessertes Verfahren zu schaffen. Die Aufgabe wird mit dem Verfahren mit den Merkmalen des Anspruchs 1 und die Messeinrichtung mit den Merkmalen des Anspruchs 11 und 12 gelöst. Vorteilhafte Ausführungsformen sind in den Unteransprüchen beschrieben. Es wird vorgeschlagen, dass zum Bestimmen der Lage eines Objektes kumulierte Drehwinkelinkremente ^^ ^^ für die Messzeitpunkte ^^ ^^ iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen ^^ ^^ die an aufeinanderfolgenden Messzeitpunkten ^^ ^^ mit k = 1, …, n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement ^^ ^^―1 und der transzendenten Funktion ^^(Ϛ) = Ϛ 1 2(1 ― Ϛ ^^ ^^ ^^ Ϛ) mit = ^^―1 ^^―1 als Funktion des ^^ /2 ⋅ ^^ /2 kumulierten Drehwinkelinkrementes ^^ ^^―1 zum jeweils vorhergehenden Messzeitpunkt ^^ berechnet werden. Hierzu erfolgt ein Approximieren der transzendenten Funktion durch Padé-Approximationen. Die Drehwinkelinkremente und Drehwinkel enthalten Winkelinformationen für die erforderlichen Raumrichtungen und können Vektoren mit den Vektorkomponenten der drei Einheiten i, j, k des oben beschriebenen Einheitsquaternions bzw. der drei Raumrichtungen x, y, z haben. Auch wenn die Drehwinkelinkremente und Drehwinkel in der vorliegenden Anmeldung und den Ansprüchen zur Vereinfachung nicht explizit als Vektoren bezeichnet sind, ist die Ausführungsform als Vektor oder Darstellungen durch schiefsymmetrische Matrizen davon umfasst. Das Verfahren hat den Vorteil, dass die Berechnung der iterativen Drehwinkel- inkremente direkt ohne Umwege über Quaternionen erfolgt. Dabei ist nur eine transzendente Funktion mit als Funktion des kumulierten Drehwinkelinkrementes ^^ zum jeweils vorhergehenden Messzeitpunkt ^^ ^^―1 zu berechnen. Dies erfolgt durch Padé-Approximationen der transzendenten Funktion. Auf diese Weise kann die einzige zu berechnende transzendente Funktion mit hoher Präzision durch eine rationale Funktion approximiert und somit effizient und mit vorab bekannter Anzahl an Rechenoperationen berechnet werden. Padé-Approximationen selbst sind bekannte Werkzeuge für numerische Approxi- mationen. Die Approximationen der transzendenten Funktion 1 ^^(Ϛ) =Ϛ2 (1 ― ^^ ^^ ^^ Ϛ) mit Ϛ als Funktion des kumulierten Drehwinkelinkrementes ^^―1 ^^ durch Approximationen nutzt die Tatsache aus, dass die Funktion ^^(Ϛ) spezielle Eigenschaften besitzt, durch die ihre Padé-Approximationen erst so genau und damit effizient werden. Damit lässt sich auch der maximale Approximationsfehler abschätzen. Zudem werden nur gerade Potenzen des Argumentes benötigt, so dass die Berechnung der Quadratwurzel der transzendenten Funktion sind nur für den Bereich 0 ≤ notwendig, da sich die Funktion für den doppelten Winkel daraus Damit lässt sich auch direkt die Singularität der Funktion bei = ^^ vermeiden, wo die Approximation zusammenbricht. Zusätzlich zu relativen kumulierten Winkelinkrementen ^^ ^^, deren Norm zur Gewährleistung einer vorgegebenen Präzision begrenzt sein muss, können auch absolute Winkelinkremente und somit insbesondere auch der die Lage des Objektes (Körpers) parametrisierende Drehwinkel mit vorgegebener Präzision berechnet werden. Die einzige dabei zu berechnende transzendente Funktion ^^(Ϛ) wird dennoch mit hoher Präzision durch eine rationale Funktion approximiert und kann somit effizient und mit vorab bekannter Anzahl an Rechenoperationen berechnet werden. Zudem bietet sich die Möglichkeit, die trigonometrischen Funktionen (Sinus und Cosinus) der kumulierten Drehwinkelinkremente ^^ ^^ und daher auch die zugehörige Rotationsmatrix ohne weitere Approximation oder Berechnung von transzendenten Funktionen aus dem Ergebnis für die transzendente Funktion ^^(Ϛ) zu berechnen. Damit lassen sich auch mit reduzierter Prozessorleistung sehr präzise und schnell trigonometrische Funktionen und Drehungen automatisiert rechnergestützt allgemein aus Winkeln ^^ ^^ aus dem Ergebnis der transzendenten Funktion ^^(Ϛ) bestimmen. Damit werden folgende technische Probleme gelöst: 1) Es können statt relativer Winkelinkremente auch absolute Winkelinkremente an einen Bordcomputer weitergegeben werden. Dadurch kann bei einem temporären Ausfall der Kommunikation zwischen IMU und Bordcomputer und somit dem Verlust auch nur einzelner Winkelinkremente gewährleistet werden, dass sich dieser Ausfall / Fehler auch ohne andere Informationsquellen (Sensoren) kompensieren lässt. Damit ist eine Lagebestimmung aus den Messdaten der Gyroskope direkt innerhalb der IMU, unabhängig vom Bordomputer möglich. 2) Die Notwendigkeit, Wertetabellen (eng. lookup tables) und/oder zeitverzögernde lterationsalgorithmen (mit ggf. einer Genauigkeitsüberprüfung als Abbruchbedingung und damit einer nicht festgelegen Anzahl an Rechenoperationen), für die Berechnung der transzendenten Funktionen zu verwenden, entfällt. Stattdessen kann eine direkte Berechnung der erforderlichen transzendenten Funktion mit vorab festgelegter Anzahl von Rechenoperationen bei garantierter Genauigkeit erfolgen. 3) Mit Hilfe der Rotationsmatrix können aus den durch die Beschleunigungssensoren der IMU gemessenen Beschleunigungen / relativen Geschwindigkeitsinkremente die translatorische Geschwindigkeit und die Position berechnet werden. In Kombination mit der Lage und mit einem Modell für die Gravitation erscheint somit eine vollständige Trägheitsnavigation (eng. inertial navigation) innerhalb der IMU möglich, unabhängig vom Bordcomputer. 4) trigonometrische Funktionen und Rotationsmatrizen können auch in anderen Anwendungsgebieten wie beispielsweise Computergrafik / Computer Vision, CAD bei gegebenen Drehwinkeln effizient und präzise berechnet werden. Die Berechnung der kumulierten Drehwinkelinkremente ^^ ^^ für einen jeweiligen Messzeitpunkt ^^ ^^ kann als Funktion des in mindestens zweiter Ordnung im zweiten Argument approximierten Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) erfolgen. Dabei werden iterativ kumulierte Drehwinkel- inkremente ^^ ^^ für die Messzeitpunkte ^^ ^^ mit k = 1,…, n jeweils aus der Funktion des approximierten Exponenten der Baker- Formel mit jeweils einem zum Messzeitpunkt ^^ gemessenen relativen Drehwinkelinkrement ^^ ^^ und dem für den ^^ vorhergehenden Zeitpunkt ^^ ^^―1 bestimmten kumulierten Drehwinkelinkrement ^^ ^^―1 oder einem vorgegebenen Start-Drehwinkelinkrement ^^ 0 als Eingangsgrößen bestimmt. Damit ist eine präzise Summation / Kumulation der Winkelinkremente auch im Fall einer zeitlich nicht konstanten Drehachse auf sehr effiziente Weise mit reduziertem Rechenaufwand möglich. Liegt zum Abtastzeitpunkt ^^ ^^ ein neues (gemessenes) relatives Winkelinkrement vor, welches abkürzend mit ^^ ^^ bezeichnet wird, so wird dieses durch ^ ^^ zu dem bestehenden kumulierten Winkelinkrement hinzugefügt und ergibt dann das neue kumulierte Winkelinkrement ^^ ^^. Dieser Ablauf wird nachfolgend in abkürzender Schreibweise als Algorithmus 1 dargestellt, wobei das Zwischenergebnis mit ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ bezeichnet ist, welches vom Zeitpunkt ^^ ^^―1 ausgehend, eine ganze Anzahl von k relativen Winkelinkrementen ^^ ^^ berücksichtigt. Das gespeicherte Winkelinkrement ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ wird nach einer Anzahl von n Updates wieder auf Null gesetzt. Damit wird ein kumuliertes, aber immer noch relatives (also auf den jeweiligen Startzeitpunkt ^^ der Kumulation bezogenes) Winkelinkrement ^^ ^ ^ ^ ^ , ^ ^^ erzeugt, ^ ^^ welches die aus der nacheinander ausgeführten Anzahl von n relativen Drehungen mit Winkeln ∆ ^^ ^^, k = 1, ...,n resultierende Gesamtdrehung parametrisiert. Algorithmus 1 zur iterativen Berechnung von ^ des Intervalls [ ^^ ^^―1 = ^^ ^^ ^^ ^ ∆ ^^ ] ^^ 1: procedure ^^ ^ ^ ^ ^ ^^ ^^ (∆ ^^ 1 , …, ∆ ^^ ^^ ) 2: Setze ^^ ^ ^ ^ ^ , ^ 0 ^ ^^ = 0 3: for k = 1,…,n do 4: Berechne ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ aus ^^ ^ ^^ ^, ^ ^ ^^ ^― ^ 1 und ^^ ^ ^^ gemäß Gleichung (1.2) Die Berechnung der Funktion ^^ in Zeile 4 des Algorithmus 1 erfordert hierbei insbesondere die lediglich approximative Berechnung nur einer einzigen transzendentalen Funktion (1.3) die mit dem Wert des letzten kumulierten Winkelinkrementes, also ^^ ^ mit auszuwerten ist (vgl. die expliziten (3.5a) weiter unten). Hierbei ist die Funktion ^^(Ϛ) gerade, so dass nur die quadrierte Norm (θ ^^―1)2  benötigt wird und damit keine Berechnung der Quadratwurzel erforderlich ist. Darüber hinaus basiert der Algorithmus auf einer einfacheren Iteration und sollte die Analyse der Fortpflanzung von Messunsicherheiten und Messrauschen gegenüber dem zum Stand der Technik einleitend beschriebenen quaternionenbasierten Algorithmus vereinfachen. Diese Eigenschaften erlauben eine direkte Implementierung auf der Hardware (zum Beispiel einem Field Programmable Gate Array (FPGA)) der IMU. Der Algorithmus kann damit als Teil eines Sensorpaketes und der Signalverarbeitung der Sensordaten für die Ausgabe an der Benutzerschnittstelle implementiert werden. Er berechnet Winkelinkremente mit reduzierter Frequenz, die im Bordcomputer dann für die Berechnung der Lage verwendet werden können. Die Padé-Approximation erlaubt eine besonders effiziente Berechnung der transzendenten Funktion im Zusammenhang mit dem iterativen Algorithmus 1. Dabei können die Drehwinkelinkremente ^^ ^^ an aufeinanderfolgenden Messzeitpunkten ^^ ^ mit k = 1, …, n gemessen (bzw. aus gemessenen ^^ ^ ^^ ^^ berechnet) und zu einem kumulierten Drehwinkel ^^ ^ für einen die Messzeitpunkte ^^ ^^ umfassenden Kumulations-Zeitbereich I zusammengefasst werden, indem jeweils ein kumulierter Drehwinkel ^^ ^ für einen Kumulations- Zeitbereich ^^ als Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) berechnet wird. Dabei werden iterativ kumulierte Drehwinkelinkremente ∆ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ für die Zeitpunkte ^^ ^^ mit k = 1,…, n jeweils aus der Funktion der Exponenten mit jeweils einem Drehwinkelinkrement ∆ ^^ ^^ des Messzeitpunktes ^^ ^^ und dem für das vorhergehende Zeitintervall [ ^^ ^^―1, ^^ ^^―1+ ^^] bestimmten kumulierten Drehwinkelinkrement ∆ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^ ― ^ 1 oder einem vorgegebenen Ausgangs- Drehwinkelinkrement ∆ ^^ ^ ^ ^ ^ , ^ 0 ^ ^^ als Eingangsgrößen bestimmt. Der kumulierte Drehwinkel ∆ ^^ ^ ^ ^ ^ ^^ ^^ für den Zeitbereich ^^ ergibt sich als Funktion des letzten kumulierten Drehwinkelinkrements ∆ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^, das für den letzten Zeitpunkt ^^ ^ ^ des Intervalls des Zeitbereichs I bestimmt wurde. Es kann diesem kumulierten Drehwinkelinkrement ∆ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ entsprechen. Für die Iteration kann das Ausgangs-Drehwinkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^ = ^ 0 für k = 0 bei einer Initialisierung mit dem Wert Null vorgegeben werden. Zur Bestimmung des kumulierten Drehwinkels ^^ ^^ ^ ^ ^ ^ ^^ ^^ für den Zeitbereich [ ^^I-1, ^^I] mit den gemessenen Drehwinkelinkrementen ^^ ^^ ^^ für die jeweiligen Zeitpunkte mit k = 1,…, n als ^^ ^^ Eingangsgrößen wird vor dem ersten Iterationsschritt die durchgeführt. Anschließend werden die kumulierten Drehwinkelinkremente ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ als Zwischenergebnisse iterativ mit k = 1 bis n aus dem jeweils für den vorhergehenden Iterationsschritt k-1 bestimmten oder vorgegebenen kumulierten Drehwinkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^ ― ^ 1 und dem für den Zeitpunkt ^^ ^^, der zum jeweiligen Iterationsschritt k gehörig ist, gemessenen Drehwinkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ als Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) berechnet. Dieser Exponent muss hierzu nur in seinem zweiten Argument durch sein Taylorpolynom approximiert werden, da der erste Exponent exakt ist. In jedem Iterationsschritt ist die Auswertung nur einer Gleichung, die später als Gleichung (3.5) näher erläutert wird, erforderlich. Dies hat die folgenden Konsequenzen: – Es ist nur eine einzige transzendente Funktion, ^^(Ϛ) = ^ zu berechnen, die effizient durch die Padé-Approximation approximiert – Die Iteration ist vereinfacht, da (bis auf die in jedem Fall notwendige Initialisierung) keine Ausführung eines zusätzlichen Rechenschrittes bei jedem n-ten Zeitschritt erforderlich ist. Dies ist insbesondere für die Implementierung auf stark limitierter Hardware (wie beispielsweise einem Field Programmable Gate Array (FPGA)) ein wichtiger Vorteil. – Die Analyse der Fortpflanzung von Messunsicherheiten und Messrauschen wird erleichtert. – Die Padé-Approximation der einzigen transzendente Funktion, ^^ θ ^^―1 2 ist genauer als die Approximation der Gleichung (1.3) durch ihr Taylorpolynom gleicher Ordnung. Für große kumulierte Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ und damit insbesondere für eine größere Anzahl n wird eine größere Genauigkeit erreicht. – Durch eine getrennte Wahl der beiden Ordnungen für die Approximation lässt sich das Verfahren recht einfach an die gegebenen Anforderungen (bspw. die maximale Drehrate, Sensorfrequenz, Teiler zur Frequenz, Anzahl der zu summierenden Winkelinkremente) für den Benutzer anpassen. Im Fall einer zeitlich veränderlichen Drehachse ergeben sich nun folgende Änderungen (im vorgeschlagenen Verfahren berücksichtigt): 1. Das von den Sensoren gemessene Integral Δφk der Drehrate ω ergibt nicht direkt das Winkelinkrement ^ das direkt Drehachse und den Drehwinkel der ^ im Abtastintervall erfolgten Drehung angibt. Stattdessen ist ^^ aus den ^^ Integralen der aktuellen und vergangenen Drehraten zu berechnen. 2. Die direkte Summation zeitlich aufeinanderfolgender mit den Zeitintervallen ^] assoziierter Winkelinkremente ^^ ^ , ∆ ^^ ergibt nicht das zur Drehung im Zeitintervall [ ^^meas, ^^―2, ^^meas, ^^] gehörende Winkelinkrement. Anstelle der einfachen Summation der Winkelinkremente ist eine kompliziertere mathematische Operation auszuführen. Das Verfahren beruht auf den bei zeitlich veränderlichen Drehachsen erforderlichen komplexen mathematischen Operationen anstelle einer einfachen Summation der Winkelinkremente. Die Ursache hierfür ist, dass Drehungen um nichtparallele Drehachsen miteinander nicht vertauschen. Vielmehr bilden die Drehungen im (dreidimensionalen) Raum eine Gruppe, die Drehgruppe SO(3), welche eine sogenannte Lie-Gruppe ist. Das mathematische Gebiet der Lie-Gruppen bildet den theoretischen Rahmen, aus welchem sich die Darstellungen der Gruppenelemente in Form von Rotationsmatrizen oder in Form der Einheitsquaternionen ergeben. Die Rotationsmatrizen sind spezielle (ihre Determinante ist 1) orthogonale (ihr Transponiertes ist ihr eigenes Inverses) Matrizen in drei Dimensionen. Der Name der Gruppe, SO(3), leitet sich aus den Begriffen Spezielle-Orthogonale Rotations- matrizen in 3 Dimensionen ab. Einheitsquaternionen sind die Erweiterung komplexer Zahlen auf drei unabhängige imaginäre Einheiten, deren aus Real- und Imaginärteil berechnete Norm 1 ist und mittels eines nichtkommutativen, normerhaltenden Produktes multipliziert werden. Eine Rotation kann durch ein Einheitsquaternion dargestellt werden, wobei die Darstellung nicht eindeutig ist, sondern sowohl q als auch -q zur Darstellung einer konkreten Drehung verwendet werden können. In Abhängigkeit der gewählten Darstellung (in Form von Vektoren, Einheitsquaternionen, …) folgen dann auch die oben beschriebenen Operationen für die Berechnung der Winkelinkremente und deren Summation. Der kumulierte Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ für einen Zeitbereich I, der eine Anzahl von n aufeinanderfolgenden Zeitpunkten ^^meas, ^^ mit k = 1, …, n beinhaltet, kann hierzu ^ ^^ ^^―1 , ^^ ^^ ^^ rekursiv (d. h. iterativ) aus einer ^ 2 2 Anzahl n von Drehwinkelinkrementen ∆ ^^ ^^, die an aufeinanderfolgenden Messzeitpunkten tmeas,k mit k = 1, …, n gemessen wurden, und dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement ^^ ^^―1 berechnet werden. Der Exponent nullter Ordnung kann dabei γ (0) ^^ ^^―1 ^^ 2 ,  =   ^^―1 ^^ 2 sein. Der Exponent erster Ordnung kann dabei ^ ^^ =   ^^ 2 ^^   +    ^ ^^ ^ ^^ ^^ wobei N 0 ein vorgegebener ganzzahliger Wert a m von Wert N0 abhängig vorgegebene der Padé-Approximation sind, und sich das für den Zeitpunkt ^^ ^^ kumulierte Drehwinkelinkrement ^ 2^ ^^ = ^^ ^^ ^^―1 ^^ 2 , jeweils aus der Summe der Entwicklungskoeffizienten γ (mα) ^^ 2 ,  des Exponenten der Baker-Campbell- Hausdorff-Formel in der nullten, ersten und mindestens zweiten Ordnung mα = 0, 1 und 2 ergibt. ^^ Es wird für den Messzeitpunkt ^^meas, ^^ ein kumuliertes Drehwinkelinkrement ^^―1 ^^ jeweils aus der Summe der Entwicklungskoeffizienten (mα) der Baker-Campbell-Hausdorff-Formel in der nullten, ersten und mindestens zweiten Ordnung mα = 0, 1 und 2 bestimmt. Die Berechnung beruht auf nur einer einzigen transzendenten Funktion, was den Zeit- und Rechenaufwand reduziert. ^^ ^^ ^^―1 ^^ ^^ ^^ kann die Summe der Argumente der Approximationen in der nullten, ersten, zweiten und dritten Ordnung mα = 0, 1, 2 und 3 umfassen. Durch die Berücksichtigung auch eines Entwicklungskoeffizienten der höheren dritten Ordnung kann die Genauigkeit weiter gesteigert werden. Auch hierbei ist nur die eine einzige transzendente Funktion enthalten, die durch die Padé- Approximation effizient und genau berechnet werden kann. Es kann weiterhin die (automatisierte, rechnergestützte) Berechnung einer Rotationsmatrix der aktiven Rotation zwischen einem ersten Koordinatensystem F des Objektes und einem nach Anwendung einer durch den kumulierten Drehwinkel ^^ ^^ beschriebenen Drehung erhaltenen zweiten Koordinatensystem G mit den Approximation von trigonometrischen Funktionen enthaltenden Ergebnissen der Padé-Approximationen und Berechnung der translatorischen Geschwindigkeit und/oder Position des Objektes aus den ebenfalls gemessenen Beschleunigungen oder relativen Geschwindigkeitsinkrementen unter Verwendung eben dieser Rotationsmatrix erfolgen. Über die Verarbeitung von Sensordaten einer IMU und der effizienten Berechnung der Lage und Position von Objekten hinaus kann eine effiziente Berechnung trigonometrischer Funktionen von gegebenen Winkeln mit reduzierter Rechenleistung und hoher Genauigkeit automatisiert rechnergestützt realisiert werden. Dies kann in vielfältige Anwendungsmöglichkeiten in verschiedenen Bereichen eingesetzt werden, z. B. in Computer Vision. Aus einem gegebenen Winkel ^^ ^^ können durch die Padé-Approximation der Funktion ^^(Ϛ) effizient und insbesondere ohne Approximation irgendwelcher weiterer transzendenter Funktionen trigonometrische Funktionen, wie beispielsweise die Funktionen ^^ ^^ 2 ^^ Ϛ2Ϛ und ^^ ^^ ^^2Ϛ, erhalten werden. Die Erfindung wird nachfolgend anhand der beigefügten Zeichnung mit einem Ausführungsbeispiel näher erläutert. Es zeigt: Figur 1 - Zeitverlauf der ^^ = gemessenen bzw. des aus den Messungen zu ^^ den Messzeitpunkten ^^meas, ^^ ermittelten Inputs ^^ ^^ ^^, dem k-ten Iterationsschritt zur Bestimmung des k-ten kumulierten Winkelinkrements Figur 2 - zur Lage eines Objektes; Figur 3 - Flussdiagramm des Verfahrens zur Bestimmung eines kumulierten Drehwinkelinkrementes aus gemessenen Drehwinkelinkrementen (der Einfachheit halber für konstante Abtastzeitschritte Δ ^^ ^^ = Δ ^^); Figur 4 - Diagramme der relativen Fehler bei der Verwendung der Padé- Approximation zur Berechnung der transzendenten Funktion ^^(Ϛ); Figur 5 - Diagramme der relativen Fehler bei direkter Verwendung der Padé- Approximation von ^^(Ϛ) und eines regularisierten Ausdruckes zur Berechnung der transzendenten Funktion ^^ Figur 6 - Relativer Fehler der Approximation von ^^ ^^ ^^ und ^^ ^^ ^^ ^^0( ^^), die auf einer Padé-Approximation von ^^(Ϛ) basieren; Figur 7 - Relativer Fehler der Approximation durch die Taylor-Polynome mit den Graden 2n+1, 2n von ^^ ^^ ^^( ^^) und ^^ ^^ ^^( ^^); Figur 8 - Relativer Fehler der Approximationen von ^^ ^^ ^^ θ und ^^ ^^ ^^ θdurch die Taylor-Polynome mit den Graden 2n+1, 2n von ^^ ^^ ^^ und verwenden in den Ausdrücken für ^^ ^^ ^^ θ und ^^ ^^ ^^ θ der halben Winkel; Figur 9 - Weiteres Flussdiagram des Verfahrens. Figur 1 zeigt einen Zeitverlauf der ^^ = ^^ gemessenen bzw. des aus den Messungen zu den Messzeitpunkten ^^ ermittelten Messwerte für die Dreh- ^ winkelinkremente, d. h. den Inputs ^^ ^^ dem k-ten Iterationsschritt zur Bestimmung des k-ten kumulierten Winkelinkrements ^^ ^^ ^ o ^, u ^^ t und dem kumulierten Drehwinkel, d. h. dem Messergebnis als Output ^^ ^^ ^ o ^, u ^^ t zum Zeitpunkt ^^ ^. Damit wird die ^ Zeitabfolge mit den Messzeitpunkten der Aufsummierung für unterschiedliche Zeitintervalle k = 1, …., n der Zusammenfassung zum niederfrequenten (Low Frequency) Messergebnis deutlich. Zur Vereinfachung der Darstellung wurden zeitlich konstante Abtastzeitschritte Δ ^^ ^^ = Δ ^^ angenommen. Die Fettschreibung der Größen bezeichnet Vektoren und die Normalschrift der Größen Skalare, wie z. B. die Norm des Drehwinkelvektors. Figur 2 zeigt eine Skizze einer Messeinrichtung 1 zur Bestimmung der Lage eines Objektes 2. Die Messeinrichtung hat mindestens ein Inertiale Messeinheit IMU, die zur Messung der Drehwinkelinkremente ^^ ^^ ^^ aufeinanderfolgender Zeitpunkte ^^ ^^ für die Erfassung rotierender (kreisender) Bewegungen der Inertialen Messeinheit IMU und des damit verbundenen Objektes 2 in den drei zueinander orthogonal stehenden Raumachsen (X- bzw. Y- bzw. Z-Achse) des kartesischen Koordinatensystems eingerichtet ist. Diese gemessenen Drehwinkelinkremente ^^ ^^ ^^ werden als Eingangswerte einer Datenverarbeitungseinheit zugeführt, die zur Berechnung eines kumulierten Drehwinkels ^^ ^^ ^ o ^ ut = ∆ ^^ ^ o ^, u ^^ t für einen Kumulations-Zeitbereich I, der aufeinanderfolgende Zeitpunkte ^^ mit i = 1,…, n umfasst, eingerichtet. Außerdem ist ^^ es optional denkbar, für die der Geschwindigkeit und ggf. Position die Beschleunigungen oder relativen Geschwindigkeitsinkremente als Komponenten entlang der drei Raumachsen zu messen und der Datenverarbeitungseinheit zuzuführen. Die Datenverarbeitungseinheit 3 kann beispielsweise ein Mikroprozessor oder Mikrocontroller sein, der ein Computerprogramm mit Programmbefehlen ausführt, die den Mikroprozessor oder Mikrocontroller veranlassen, die Verfahrensschritte des anspruchsgemäßen Verfahrens auszuführen. Die Datenverarbeitungseinheit 3 kann aber auch als festverdrahtete Hardware-Logik, wie beispielsweise ein Field- Programmable-Gate-Array o. ä. ausgeführt sein. Mit dem Verfahren wird eine direkte Berechnung des kumulierten Drehwinkels ^^ ^^ iterativ aus jeweils dem im letzten Schritt gespeicherten und dem aktuell gemessenen Drehwinkelinkrement ^^ ^^ ^^ durchgeführt. Die dabei ebenfalls in diesem Verfahren erhaltene Rotationsmatrix kann zur Berechnung der Geschwindigkeit und ggf. Position aus den ebenfalls gemessenen Beschleunigungen oder Geschwindigkeitsinkrementen verwendet werden. Grundlage für die Summation der gemessenen Drehwinkelinkremente bildet die aus der Theorie der Lie-Gruppen stammende Baker-Campbell-Hausdorff (BCH) Formel, die ein Vertauschungsgesetz für bestimmte lineare Operatoren angibt. Sie liefert für die Multiplikation zweier durch die Exponentialabbildung erhaltenen Elemente eα, eβ einer Lie-Gruppe den Exponenten γ(α, β) des resultierenden Elementes e γ(α, β), welches also die Relation eγ(α, β) = eα eβ (3.1) erfüllt. Das Ergebnis γ(α, β) wird durch die Baker-Campbell-Hausdorff Formel als formale Reihenentwicklung in der Gesamtordnung m von ^^ und ^^ angegeben und hat daher schematisch folgende Form: ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ mit dem Koeffizienten ^^( ^^ ^^, ^^ ^^), wobei die tatsächliche Form wegen der allgemeinen Nichtvertauschbarkeit der Gruppenelemente ^^ ^^ ≠ ^^ ^^ wesentlich komplizierter ist. Die allgemeine Baker-Campbell-Hausdorff (BCH) Formel ist als formale Reihe noch nicht in einer Form, die für die Summation zeitlich aufeinanderfolgender mit den Zeitintervallen [ti−3+k, ti-2+k], [ti-2+k, ti-1+k] assoziierter Winkelinkremente ∆ ^^ ^^―1 , ∆ ^^ ^^ geeignet ist, um das tatsächliche zur Drehung im Zeitintervall [ti−3+k, ti-1+k] gehörende Drehwinkelinkrement, d. h. den Drehwinkel ^^ ^^ ^ o ^ ut für den Zeitbereich I (in diesem Fall bestehend aus zwei Unterintervallen), zu erhalten. Die BCH-Formel muss noch auf den Spezialfall der Lie-Gruppe SO(3) spezialisiert werden. Zudem muss eine geeignete Auswertung / Berechnung der Reihenentwicklung erfolgen. Nach der Spezialisierung auf SO(3) wird eine geeignete Auswertung des dann immer noch als Reihe vorliegenden Ausdrucks durchgeführt. Die Spezialisierung auf die Lie-Gruppe SO(3) liefert zunächst nur die benötigte Relation für die Summation zweier Winkelelemente, für den zuvor beschriebenen Fall ^ ^^ wobei γ nach wie vor als Reihenentwicklung in der Gesamtordnung m wie in Gleichung (3.2) gegeben ist und ^^ ^^ das neu kumulierte Drehwinkelinkrement ist, das aus dem neu gemessenen relativen Drehwinkelinkrement ∆ ^^ ^^ und dem vorhergehend berechneten kumulierten Drehwinkelinkrement berechnet wird. ^^ Die Auswertung der Reihenentwicklung gemäß Gleichung (3.2) der Funktion γ in der Gleichung (3.3) erfolgt unter der Annahme, dass nur ihr zweites Argument klein ist und daher nur in diesem Argument eine Approximation durch das Taylorpolynom vorgenommen werden kann. Dieses ist erforderlich, denn bei iterativer Anwendung der Relation (3.3) auf insgesamt n Elemente gemäß der Vorschrift: ^ ^ ^^ beinhaltet das erste Argument von γ das (iterative) kumulative Drehwinkelinkrement ^ aus den vorangegangenen k–1 Schritten. Dieses ist also für große k ≤ n Weise klein. Daher kann eine Approximation der Taylorreihe im ersten Argument durch ihr Taylorpolynom zu einem inakzeptabel großen numerischen Fehler führen. Die BCH-Formel lässt sich nun für die Lie-Gruppe SO(3) wie folgt approximieren, um die Exponenten ^^ ^^ der Ordnungen mβ = 0, 1, 2 zu erhalten: ^ ^^ ^ ^^ mit der transzendenten Funktion Die einzige, in dem Gleichungssystem enthaltene und durch Approximation zu berechnende Funktion ist die transzendenten Funktion ^^(Ϛ). Die Herleitung der Approximation mindestens zweiter Ordnung mβ der Baker- Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) wird später noch im Detail beschrieben. Es reicht aus, wenn lediglich im zweiten Argument ^^2 = ^^ ^^ eine Approximation mindestens in den Ordnungen mβ = 0, 1, 2 vorgenommen wird. Die Ausdrücke sind hingegen im ersten Argument ^^ ^ exakt. Es tritt nur eine einzige transzendente Funktion B θ 21   = B θ ^^―1 2 auf, die mit Hilfe der Padé-Approximation effizient und genau approximiert werden kann. ^ und den von bei gegebenen maximal zu erwartenden Werten für die Winkelinkremente und deren Anzahl n ab, die in einem einzigen Winkelinkrement ^^ ^^ ^ ^ ^ ^ ^^ ^^ zusammengefasst werden sollen. Bei Verwendung der Gleichungen (3.5) kann die erforderliche Ordnung bei gleicher numerischer Genauigkeit deutlich niedriger gewählt werden als bei Verwendung der Gleichungen (2.3), (2.5) und (2.6). Damit kann die Komplexität reduziert und der Rechenaufwand, d. h. Berechnungszeit und Hardwareaufwand, reduziert werden. Die iterative Berechnung des kumulierten Drehwinkels ^^ ^^ ^ ^^^ ^^ ^^ des Intervalls [ ^^I-1, ^^I] mit Hilfe der Approximation mindestens zweiter Ordnung der Baker- Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) kann auf Grundlage der Gleichungen (3.4) und (3.5a), (3.5) durch die Prozedur 1 wie folgt durchgeführt werden: Prozedur ^^ ^^ ^^ Setze Für ^^ iterativ: Berechne ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ aus dem (iterativen) kumulierten Drehwinkelinkrement ^^ ^^ ^ ^^ ^, ^ ^ ^^ ^― ^ 1 und dem gemessenen Drehwinkelinkrement ^^ ^^ ^^ gemäß Gleichung (3.4) mit der expliziten Form (3.5a), (3.5) Setze den resultierenden kumulierten Drehwinkel ^^ ^^ ^ ^ ^ ^ ^^ ^^ gleich dem letzten (iterativen) kumulierten Drehwinkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^. Figur 2 zeigt ein zugehöriges Flussdiagramm dieses Verfahrens bzw. der obigen Prozedur zur Bestimmung eines kumulierten Drehwinkelinkrementes aus mit einem Drehwinkelsensor für aufeinanderfolgende Zeitintervalle gemessenen Drehwinkel- inkrementen. Berechnung des expliziten Ausdrucks für die BCH Formel für die Lie-Gruppe SO(3) Im Folgenden wird ein Ausdruck für die BCH Formel für die Lie-Gruppe SO(3) hergeleitet, die Rotationen im dreidimensionalen Euklidischen Raum beschreibt. Strategie ist hierbei, die Entwicklung der Formel in Ordnungen nur ihres zweiten ^^2 = ^^ ^^ ^ ^^^ ^^ ^^ durchzuführen und die Abhängigkeit vom ersten Argument ^^ 21 zu resummieren, so dass die Formel bezüglich ihres ersten Argumentes ^ 2^1 exakt ist. Das Ergebnis stellt also eine Entwicklung der Formel in Ordnungen ihres zweiten Argumentes bei gleichzeitiger Resummation aller Abhängigkeiten von ihrem ersten Argument dar. M. Müger. „Note on the theorem of Baker-Campbell-Hausdorff-Dynkin“, https://www.math.ru.nl/~mueger/PDF/BCHD.pdf gibt einen Überblick über einige Grundlagen, Fakten und Beweise zur BCH-Formel. Dort wird im Abschnitt 8 und insbesondere in Bemerkung 8.2 auf Seite 24 darauf hingewiesen, dass auf Grund- lage einer speziellen Stufung der Algebra, basierend auf der Wortlänge im zweiten Argument anstelle der Gesamtwortlänge in beiden Argumenten, eine Entwicklung der BCH Formel im zweiten Argument erfolgen kann. Der entsprechende Ausdruck ist bis zu ersten Ordnung im zweiten Argument angegeben, ohne allerdings eine Spezialisierung auf die Lie-Gruppe SO(3) vorzunehmen und die die Bernoulli-Zahlen beinhaltende Serie explizit zu resummieren. 1. Berechnung der nullten und ersten Ordnung, ^^ ^^ = 0,1 Wird die in erster Ordnung einfachen Schritte der Spezialisierung und der Resummation unter Verwendung der Erzeugendenfunktion der Bernoullizahlen aus, welche durch ^ ^ gegeben ist, ergibt sich für den in M. Müger „Note on the theorem of Baker- Campbell-Hausdorff-Dynkin“. Als Reihe gegebenen Ausdruck der BCH-Formel die folgende explizite Darstellung in resummierter Form: enthält nur gerade Ordnungen von so dass die Berechnung der Gleichung (A.2) keine Berechnung der in der Norm ^^ des Vektors enthaltenen Quadratwurzel erfordert. Dieses ist bei numerischen Vorteil. 2. Grundlagen zum Formalismus und Zusammenhang zwischen der Differentialgleichung der Winkelinkremente und der BCH Formel Die zur Herleitung der zweiten Ordnung erforderlichen Grundlagen und Schritte werden im Folgenden dargestellt. Ausgangspunkt bildet die Exponentialabbildung, die auch für Argumente α definiert ist, die Elemente einer Lie-Algebra, also des Tangentialraums einer Lie-Gruppe am Element der Identität sind, und die durch Matrizen dargestellt werden können. Die Exponentialabbildung ist über ihre Taylorreihe definiert ^ ^ wobei α ein Element einer Darstellung einer Lie-Algebra ist. Dieses Element α kann z. B durch die reinen Quaternionen im Fall der Rotationen gegeben sein. Potenzen in der obigen Gleichung sind so zu verstehen, dass sie mit dem entsprechenden Produkt durchzuführen sind, was für die Darstellung definiert ist, also z. B. dem Matrixprodukt im Fall einer Matrixdarstellung oder dem Quaternionenprodukt im Fall einer Darstellung durch Quaternionen. Da verschiedene Elemente α, β einer Lie-Algebra in der Regel nicht miteinander vertauschen, also ^^ ^^ ≠ ^^ ^^ gilt, ist die Ableitung der Exponentialabbildung nicht durch den bekannten Ausdruck gegeben. Vielmehr liefert die Anwendung der Ableitung auf die Reihendarstellung gemäß Gleichung (A.4) den folgenden Ausdruck: wobei der n-fache Kommutator folgendermaßen definiert ist: und die adjungierten Wirkung des Elements α auf ein anderes Element β repräsentiert. Aus der Ableitung gemäß Gleichung (A.5) folgt nun die Relation (A.7) Diese lässt sich invertieren. Dazu wird die Erzeugendenfunktion der Bernoulli-Zahlen verwendet, welche in Gleichung (A.1) gegeben ist. Das Inverse von der Relation (A.7) ist dann gegeben durch Es wird nun eine für numerische Berechnungen in der Praxis bedeutsame Eigenschaft erkennbar, nämlich, dass die Reihe in der Inverse gemäß (A.8) nur gerade Potenzen von α enthält. Deshalb hängt der dann auf die Lie-Gruppe SO(3) spezialisierte Ausdruck nur von geraden Potenzen der Norm von α ab und beinhaltet somit keine Quadratwurzeln. Zunächst wird die Relation für die Erzeugendenfunktion auf ihre geraden und ungeraden Potenzen projiziert, so dass die Gleichungen erhalten werden: Dieses Ergebnis wird verwendet, um die Inverse gemäß Gleichung (A.8) in folgende Form zu bringen: (A.11) Dieser Ausdruck (A.11) wird nun spezialisiert auf die Lie-Gruppe SO(3), bzw ihrer doppelten Überlagerung, die Lie-Gruppe SU(2), dargestellt durch die Einheitsquaternionen. Für zwei Elemente θ, ω aus ihrer Lie-Algebra ^^ ^^(2) reduzieren sich die multiplen Kommutatoren wie folgt: Durch wiederholte Anwendung dieser Relation ergibt sich für die Gleichung (A.11) dann der Ausdruck Nach Identifikation des einfachen Kommutators mit dem Kreuzprodukt und Reduktion des doppelten Kommutators / Kreuzproduktes ist das Ergebnis (A.13) in der Tat die Differentialgleichung für das Winkelinkrement θ in Abhängigkeit von der Winkelgeschwindigkeit ω, die in der Gleichung (2.1) gegeben ist. Das Ergebnis der Gleichung (A.8) und im Spezialfall der Lie-Gruppe SO(3) dessen expliziter Ausdruck durch (A.13) sind eng verknüpft mit der Baker-Campbell- Hausdorff Formel, wie im Folgenden gezeigt wird. Die Baker-Campbell-Hausdorff Formel liefert das Lie-algebra-wertige Element γ als Funktion zweier Lie-algebra- wertiger Elemente α, β, als Lösung folgender Relation ^^ ^^( ^^, ^^) = ^^ ^^ ^^ ^^ . (A.15) Um diesen Zusammenhang zu zeigen, wird der Exponent γ in Ordnungen von β entwickelt, wie dieses z. B. in M. Müger „Note on the theorem of Baker-Campbell- Hausdorff-Dynkin“. Dargestellt ist. Hierzu wird einen Parameter t ∈ R eingeführt und die Entwicklung geschrieben als wobei klarer Weise ^^ (0) = ^^( ^^, 0) = ^^ ist. Das Ergebnis für ^^ in der Entwicklung der Gleichung (A.16) kann direkt aus der Gleichung (A.11) werden, wenn α durch γ(α, t β) ersetzt wird. Dieses führt dann auf die folgende Differentialgleichung für γ: Durch Einsetzen von t = 0 in die obige Gleichung wird wegen der Gleichung (A.16) ein Ergebnis für ^^ Ein Vergleich dieses Ergebnisses mit der Gleichung (A.11) zeigt, dass die Differentialgleichung für das Lie-algebra-Element α bei gegebenem Lie-algebra- Element β = e−α ∂t eα identisch mit der Gleichung für den in β linearen Beitrag in der BCH Reihe gemäß Gleichung (A.16) ist. Hinzu kommt, dass durch Einsetzen von t = Δt die Approximation von γ durch ^^ ≅ ^^ + ^^ (1) ( ^^, ∆ ^^ ^^) + ^^(∆ ^^ ^^) 2 zu einer Lösung der Differentialgleichung (A.11) zur linearen Ordnung in Δt β wird. 3 Berechnung der zweiten Ordnung, ^^ ^^ = 2 Im Folgenden wird der Term ^^ (2) der BCH Reihe (A.16) hergeleitet. Dieser beinhaltet die Terme 2. Ordnung in β. Zunächst müssen hierfür die Ableitung der Relation (A.11) berechnet werden. Dies führt zu: Eine Gleichung für γ(2) ist analog zum Vorgehen im Fall für γ(1) zu erhalten. Hierzu kann α durch γ(α, t β) ersetzt und außerdem ^^ = 0, ^^ = 0 eingesetzt werden. Damit ergibt sich Um ein explizites Ergebnis für den Spezialfall der Lie-Gruppe SO(3) zu erhalten, wird die Reihe in (A.20) zu einem geschlossenen Ausdruck summiert. Es kann analog wie zur Herleitung von (A.13) vorgegangen werden, wobei die Gleichung (A.12) verwendet wird, um die multiplen Kommutatoren zu reduzieren. Außerdem wird direkt ^^ = ^ 2^1 , ^^ = ^ 2^2 gesetzt, um das Ergebnis in den im Haupttext verwendeten Variablen zu erhalten. Um die Reduktionsformel (A.12) nutzen zu können, kann zunächst festgestellt werden, dass sich eine allgemeine Folge mit Argument fk wie folgt in Folgen der nur geraden und nur ungeraden Anteile aufteilen lässt: (A.21) Mit dieser Identität lässt sich die innere Summe in (A.20) wie folgt aufteilen: wobei ϑ die Heavyside Stufenfunktion ist. Einsetzen dieses Ergebnisses in die Reihe in (A.20) liefert die Taylorreihe: Ein geschlossener Ausdruck für die erste der beiden Reihen ist bereits aus (A.9) bekannt. Ein geschlossener Ausdruck für die zweite Reihe kann durch Anwenden eines geeigneten Ableitungsoperators auf das bekannte Ergebnis erhalten werden: Analytische Fortsetzung der Ausdrücke durch Identifikation z  =  i θ1 führt dann für den Term der BCH Reihe, der quadratisch in θ2 ist, auf das Ergebnis das nur von einer einzigen transzendenten Funktion, B(ζ), abhängt, die bereits im Ausdruck für γ (1) erscheint und in (A.3) gegeben ist. Dieser Umstand ist für numerische Berechnungen von Vorteil, da nur eine einzige transzendente Funktion approximiert und berechnet werden muss. 4. Expliziter Ausdruck für die BCH Formel bis zur zweiten Ordnung, mβ = 0, 1, 2 Das explizite Ergebnis der BCH Formel für SO(3) in der Entwicklung (A.16), d. h. bis einschließlich der Terme quadratischer Ordnung im zweiten Argument ^ 2^2 wird im Folgenden zusammengefasst und vereinfacht. Der Term ist einfach durch das erste Argument ^^1 der Funktion gegeben. Der Term ist durch (A.13) gegeb 2 en, wenn darin ^^ →  ^^1 und ^^ →  ^^2 ersetzt oder wahlweise direkt der Ausdruck (A.2) verwendet wird. Der Term γ (2) ist in (A.25) dargestellt. Damit ergibt sich für die BCH Formel Nun werden die folgenden Identitäten verwendet: um zunächst die in der Gleichung (A.26) auftretenden multiplen Kommutatoren zu berechnen und zu vereinfachen Damit kann die Formel (A.26) schrittweise auf folgende Form gebracht werden: 5. Berechnung höherer Ordnungen, mβ ≥ 3 Das beschriebene Verfahren lässt sich direkt auch für die Berechnung höherer Ordnungen anwenden. Insbesondere die Ordnung mβ = 3 in θ2, also der Term γ (3) in der BCH Reihe (A.16) kann genutzt werden, um eine noch größere Genauigkeit zu erhalten. Die Berechnung der Exponenten γ (mα) in der nullten, ersten und mindestens zweiten Ordnung mα = 0, 1 und 2 erfordert eine approximative Berechnung lediglich einer einzigen transzendenten Funktion ^^(Ϛ) = Ϛ 1 2(1 ― ^^ ^^ ^^ Ϛ). Diese ist mit dem Wert der halben quadrierten Norm des letzten kumulierten Winkelinkrementes, also = ^ ^^ mit θ ^^―1 =   ^^ ^^―1 ⋅ ^^ ^^―1 auszuwerten. Hierbei ist die Funktion ^^(Ϛ) gerade, so dass nur die quadrierte Norm ( ^^ ^^―1 ) 2 benötigt wird und damit keine Berechnung der Quadratwurzel erforderlich ist. Für hinreichend kleine 0 ≪ 1, also kleine Norm (θ ^^ ) 2 ≪ 1 des kumulierten Winkelinkrementes ^^ ^^ kann die Funktion ^^(Ϛ) in der Gleichung (1.3) durch ihr Taylorpolynom in niedriger Ordnung approximiert werden. Dieser Fall relativ kleiner kumulierter Winkel tritt beispielsweise auf, wenn lediglich eine endliche Anzahl n relativer Winkelinkremente in ∆ ^^ ^^ kumuliert werden soll, deren Norm ebenfalls jeweils klein ist, also ^^ ^^ ^^ ≪ 1 erfüllt. Da nach n Zeitschritten der kumulierte Winkel ausgegeben und dann dessen intern gespeicherter Wert wieder auf 0 gesetzt wird, handelt es ich bei den kumulierten Winkelinkrementen ^^ ^^ ebenfalls um relative Winkelinkremente, die auf ihren jeweiligen Startzeitpunkt der Kumulation bezogen sind. Die Genauigkeit der Approximation der Funktion ^^(Ϛ) durch ihr Taylorpolynom wird somit dadurch bestimmt, wie groß n und die zu summierenden relativen Winkelinkremente ^^ mit k = 1, ...,n sind. Beispielsweise ist für n relative Winkelinkremente ^^ ^^ konstanter Norm ∆ ^^ 1 = • • • = ^^ = die Größe des ^^ ∆ ^^ Produktes n∆ ^^ ausschlaggebend für die Genauigkeit der für die Auswertung von ^^ (Ϛ) aus Gleichung (1.3) verwendeten Approximation. Im Folgenden wird die Renormierung der Winkelinkremente erläutert. Um den absoluten Drehwinkel und somit die Lage des Körpers im Raum zu erhalten, sind die ab einem Zeitpunkt, an dem die Lage bekannt war, alle relativen Winkelinkremente ^ ^^ auzuaddieren. Dieses stellt also eine Kumulation beliebig ^ vieler relativer Winkelinkremente ^^ ^^ ohne periodisches Zurücksetzen des Ergebnisses dar. Die Zahl n und daher auch das Argument Ϛ, mit welchem die Funktion ^^(Ϛ) in (1.3) auszuwerten ist, sind dann a priori nach oben nicht beschränk. Die Approximation von (1.3) durch das Taylorpolynom führt dann bei größer werdenden zu einem immer größer werdenden und letztendlich unbeschränkt und somit inakzeptablen Fehler. Dieses Problem der nach oben nicht beschränkten Norm der Winkelinkremente und somit eines unbeschränkten Argumentes Ϛ kann wie im Folgenden dargestellt gelöst werden. Dieses ist eine der Voraussetzungen dafür, die periodische, transzendente Funktion (1.3) mit nach oben begrenztem Genauigkeitsverlust durch z. B. ihr Taylorpolynom zu approximieren. Es wird darauf hingewiesen, dass eine Rotation in zwei Dimensionen durch einen einzigen Drehwinkel ϕ bestimmt ist, der auf das Intervall [-π, π] oder bei Angabe des Vorzeichens ± der Drehrichtung auf das Intervall [0, π] eingeschränkt werden kann. Das bedeutet, dass der Drehwinkel ϕ und alle Drehwinkel ϕ ^^ = ϕ ± 2π ^^, ^^ ∈ ^^ dieselbe Rotation beschreiben. Sollte also ein Drehwinkel ϕ > π angegeben sein, so ist es immer möglich, diesen durch Subtraktion geeigneter ganzzahliger Vielfache von 2π auf seinen Repräsentanten ϕ′ = ϕ ―2πκ0 , der innerhalb des Intervalls [-π, π] liegt, zu reduzieren. Wird bei gegebenem ϕ > 0 die Größe κ0 = ^^ ^^ ^^(ϕ) | 2ϕ π| gewählt, wobei [ ] die Floor-Funktion ist, dann wird mit ϕ′ = ϕ ―2πκ0 in der Tat der Repräsentant erhalten, der die Bedingung erfüllt. Bei einer Rotation in drei Dimensionen ist ein ähnliches Vorgehen möglich. Eine solche Rotation wird durch einen vektorwertigen Drehwinkel ^^ parametrisiert, dessen Norm (Länge) ^^ = ^^ ∙ ^^ den totalen (skalaren) Drehwinkel und dessen Richtung (Einheitsvektor) ^^ ^^ = ^ ^^ ^ die Drehachse festlegen. Da ^^ > 0 ist, legt ^^ ^^ auch gleich die Drehrichtung fest. Hierbei kann die Rechte-Faust- oder Korkenzieher-Regel angewendet werden, bei der der Daumen der rechten Hand in Richtung des Vektors zeigt und die Fingerspitzen der zur Faust geballten rechten Hand die Drehrichtung für Drehungen mit positivem, Winkel anzeigen, d. h. die Norm des vektorwertigen Drehwinkels. Analog zu den Rotationen in zwei Dimensionen kann daher in drei Dimensionen jede Rotation durch einen vektorwertigen Drehwinkel ^^′ = ^^′ ^^ ^^′ parametrisiert werden, dessen Norm 0 ≤ ^^′ ≤ ^^ erfüllt. Wenn ein allgemeiner Drehwinkel ^^ gegeben ist, dessen Norm ^^ nicht im Intervall [0,π] liegt, kann dieser durch einen Repräsentanten ^^′ mit minimaler Norm 0 ≤ ^^′ ≤ ^^ ersetzt werden, der aus ^^ wie folgt berechnet werden kann ^^ ′ =  (θ ― 2πκ0) ^^ ^^ (3.1) Hierbei ist ^^0 ∈ Norm und Richtungsvektor des Repräsentanten mit minimaler Norm sind dann durch ^^′ =  |θ ― 2πκ0|,      ^^θ′ = sgn(θ ― 2πκ0) ^^ ^^ (3.2) gegeben, d. h. die Bedingung 0 ≤ ^^′ ≤ ^^ ist erfüllt. Im nächsten Schritt wird der Wert bestimmt, der für κ0 zu wählen ist, um den Repräsentanten mit minimaler Norm zu erhalten. Ein allgemeiner Drehwinkel ^^ besitzt eine Norm ^^, die entweder im Intervall 2 ^^ ^^ ≤ ^^ < (2 ^^ + ^^ oder im Intervall (2 ^^ + 1) ^^ ≤ ^^ < 2( ^^ + 1) ^^ ^^, ^^ ∈ ℕ, liegt. In ersterem und letzterem Fall führt die Wahl der Konstanten κ0 = ^^ bzw. κ0 = ^^ + 1 auf folgende Ausdrücke: Diese können auch wie folgt dargestellt werden In der Tat erfüllt die Norm ^^′ nun die Bedingung 0 ≤ ^^′ ≤ ^^. Die beiden Fälle lassen sich mit Hilfe der Floor-Funktion [ ] kompakt zu Ein Vergleich der Ausdrücke aus den Gleichungen (3.5) und (3.1) ergibt, dass für einen Drehwinkel ^^, dessen Norm im Intervall κπ ≤ θ < (κ + 1)π liegt, die Konstante κ0 = κ + 21 gewählt werden muss, um dann mittels der Gleichung (3.1) den Repräsentanten ^^′ mit minimaler Norm zu erhalten. Das lässt sich unmittelbar verifizieren, wenn die Norm des Repräsentanten ^^' unter Zuhilfename der Relation ^^ ― 2 ^ 2^ = 1 ― ( 2 ―1) ^^ , welche für beliebige ^^ ∈ ℕ gültig ist, berechnet wird. Damit wird die Bedingung 0 ≤ ^^′ ≤ ^^ erfüllt. Bei dem oben detailliert beschriebenen Algorithmus für die Kumulation relativer Winkelinkremente liegt zu jedem Abtastzeitpunkt ^^ ^^ das Quadrat der Norm des gespeicherten kumulierten Winkels oder Winkelinkrementes ^^ ^^―1 vor, welches hier wieder abkürzend mit ( ^^ ^^―1 ) 2 bezeichnet wird. Auch, wenn die Norm des Winkelinkrements ^^ ^^―1 minimal ist, also die Bedingung 0 ≤ erfüllt ist, ist ^^ ≤ ^^ dieses nicht notwendiger Weise mehr für das kumulierte Drehwinkelinkrement ^^ ^^ der Fall, das mit dem auf der Gleichung (1.2) beruhenden Algorithmus und dem aktuellen relativen Winkelinkrement ^ ^^ berechnet wurde. ^ Unter der Annahme, dass ^^ ^^―1 bereits bei Bedarf durch seinen Repräsentanten mit minimaler Norm ersetzte Winkel ist, also die Bedingung 0 ^^ erfüllt ist, und ≤ ^^ dass das aktuelle relative Winkelinkrement ^^ ^^ die Norm κ ^^π ≤ Δθ ^^ < (κ ^^ + 1)π besitzt, also beispielsweise für 0 ≤ ^^ ^^ < ^^ die Konstante κ ^^ = 0 ist, für ^^ ≤ ∆ ^^ ^^ < 2 ^^ dann κ ^^ = 1 ist usw, ergibt sich für die quadrierte Norm bis einschließlich zur ^^ wobei eine Obergrenze mittels der Dreiecksungleichung für Vektoren ermittelt werden kann. Ein Vergleich dieser Abschätzung mit der Gleichung (3.4) zeigt, dass die Norm des aktuellen relativen Winkelinkrementes den Wert der Konstanten κ bestimmt, der für die Berechnung des Repräsentanten minimaler Norm erforderlich ist. Im Fall, dass alle relativen Winkelinkremente eine hinreichend kleine Norm haben, also ^^ ^^ < für alle k gilt, sind alle und damit ergeben sich in der ^^ κ ^^ = 0 Gleichung (3.5) die Fälle κ = 0, in denen das kumulierte Winkelinkrement bereits minimale Norm hat oder κ = 1. In der Praxis sollte die Einschränkung ^^ ^^ < ^^ die Regel sein, denn üblicherweise wird die Abtastzeit der Gyroskope angepasst an die zu erwartende Systemdynamik so gewählt, dass nur eine geringe Drehbewegung während eines einzigen Abtastzeitschrittes erfolgt, also ^^ ist. ^^ ≪ ^^ Im generischen Fall, wenn also individuelle κ ^^ so zu wählen sind, dass κ ^^π ≤ Δθ ^^ ^ gilt, gibt es z. B. folgende zwei Möglichkeiten für die Bestimmung von ^ κ im (k + 1)ten Zeitschritt: — Zunächst erfolgt eine Bestimmung von (3.7) Dann erfolgt eine Berechnung von κ als die natürliche Zahl, die kleiner oder gleich der Quadratwurzel von κs ist, also Mit der so erhaltenen κ Diese gibt den relativen, also auf Rest an, der nach Subtraktion der maximal möglichen ganzzahligen Vielfachen von ^^2 von dem Quadrat der Norm ( ^^ ^^)2 noch übrig bleibt. Aufgrund der Konstruktion gilt in jedem Fall 0 ≤ ^^. Mit diesem Ergebnis lässt sich dann der Repräsentant minimaler Norm durch Umformung der Gleichung (3.5) wie folgt konstruieren — Alternativ, um die (wenn auch nur approximative) Berechnung der Quadratwurzel der Konstanten κszu vermeiden, können ganzzahlige Werte κ’ über das Intervall 1 ≤ κ′ iteriert werden, entweder vorwärts startend von der unteren oder ^^ von der oberen Intervallgrenze. Dann kann mit dem aktuellen Wert κ’ jeweils berechnet und der Vorzeichenwechsel von der Größe εκ′ detektiert werden. Der größte Wert von κ’, für den εκ′ ≥ 0 gilt, ist dann der Wert für κ und εκ ≥ 0. Mit diesen Werten kann dann über die Gleichung (3.10) der Repräsentant mit minimaler Norm berechnet werden. Der Ausdruck (3.10) für die Berechnung des Repräsentant mit minimaler Norm erfordert die Berechnung der Quadratwurze ^^ 1 + ^^. Im Allgemeinen ist κπ ≤ θ ^^+1 < (κ + 1)π, woraus folgt, dass ^^ dann die Bedingung 0 ≤ ε < 2 1 κ erfüllt. Somit ist ^^ nicht notwendiger Weise klein, so dass die Quadratwurzel nicht einfach durch ihr Taylorpolynom in fester Ordnung approximiert werden kann. Jedoch wird, wie oben bereits erwähnt, in der Praxis überwiegend der Fall auftreten, dass die Norm der relativen Winkelinkremente ^^ ^^ die Bedingung ∆ ^^ ^^ erfüllt. Dieses ist ohnehin ≪ ^^ eine Bedingung, die für die Anwendbarkeit des auf dem Baker-Campbell-Hausdorff- Formel beruhenden Algorithmus gegeben sein muss und stellt damit keine zusätzliche Einschränkung dar. Dann ergibt sich aus der Gleichung (3.6), dass die Norm des kumulierten Winkels ^^ ^^ auf ^^ ^^ ≤ ^^ ^^―1 ^^ ^^ ≤ ^^ + ∆ ^^ beschränkt ist. ^^ Eine Berechnung von ^^ z. B. mittels (3.11) mit κ = κ′ = 1 ergibt dann, dass ^^ folgendermaßen beschränkt ist Wenn nun alle Zeitpunkte ^^ ^^ z. B. aus den Spezifikationen der konkreten Anwendung abgeleitet werden kann, man also ∆ ^^ so wählt, dass ^^ ≤ ∆ für alle k gilt, dann ergibt sich als Begrenzung für ^^ die Bedingung 0 ≤ ^^ ≤ 2 ∆ ^^ ^^ . Durch Erhöhung der Abtastrate können die ^^ ^^ und somit auch ∆ verkleinert werden. Eine ^^ Approximation von 1 + ^^ durch ihr Taylorpolynom wird möglich. Die maximale Ordnung der Taylorentwicklung, also der Grad des Taylorpolynoms zur Erreichung einer gewünschten Genauigkeit kann bestimmt werden. Aus der Gleichung (3.5) ergibt sich wobei das Restglied unter Verwendung der Stirlings Formel für das asymptotische Verhalten der Fakultäten wie folgt approximiert werden kann Der durch Verwendung des Taylorpolynoms vom Grad n0 verursachte Fehler in der Norm lässt sich nun wie folgt abschätzen Die Berechnung des Repräsentanten minimaler Norm kann mit folgenden Verfahrensabläufen durchgeführt werden: Algorithmus 2: Bestimmung des Repräsentanten minimaler Norm zu ^^ ^^ Es wird angenommen, dass ( ^^ ^^ ) 2 vorliegt und die Bedingung 0 ≤ ^^ ^^ ≲ ^^ + ∆ ^^ erfüllt ist. 1: procedure ^^ 2: Berechne ^^ = ( ^ ^^ ^ ―1 3: if ^^ > 0 then 4: Berechne ^^′ ^^ mittels Gleichung (3.13) (für K = 1 und dem aus der Gleichung (3.15) bestimmten n0) zu 5: else 6: Setze ^^′ ^^ = ^^ ^^ Daraus lässt sich leicht n0 bei vorgegebener Genauigkeit und vorgegebenem ∆ ^^ berechnen. Bei Überprüfung der Norm des kumulierten Winkelinkrements und ggf. Ersetzung durch den Repräsentanten minimaler Norm in jedem Zeitschritt ist im Fall relativer Winkelinkremente, die die Bedingung 0 ≤ ∆ ^^ ^^ < ^^ erfüllen, der Fall κ = 1 der einzige, in dem die Ersetzung des kumulierten Winkelinkrementes durch seinen Repräsentanten minimaler Norm durchzuführen ist. Daher lässt sich der vereinfachte Algorithmus 2 formulieren. Im Fall, dass die Norm ^^ nicht beschränkt ist, kann der allgemeinere Algorithmus ^^ 3 eingesetzt werden. Algorithmus 3: Bestimmung des Repräsentanten minimaler Norm zu ^^ ^^ Es wird angenommen, dass ( ^^ ^^ ) 2 vorliegt. 1: procedure ^^′ ^^ ( ^^ ^^ ) 2 , ^^ ^^ 2: Berechne ^^ = ( ^^ ^^)2 ^^2 3: if ^^ > 0 then 4: Berechne κ ^^ = | ^^| 5: Berechne κ = [ κ ^^] 6: Berechne ^^′ ^^zu 7: else 8: Setze ^^′ ^^ = ^^ ^^ Eine Alternative kann z. B. auch die als Algorithmus 4 gegebene Variante sein. Modifikationen der Berechnung von ^^′ ^^, z. B. in Form einer Ersetzung von ^^ durch ^^ = ^^ ― 1 und anschließender Auswertung der Taylorreihe aus der Gleichung (3.13) bis zur gewünschten Ordnung n0, können ebenfalls bei Bedarf leicht vorgenommen werden. Algorithmus 4: Bestimmung des Repräsentanten minimaler Norm zu ^^ ^^ Es wird angenommen, dass ( ^^ ^^ ) 2 vorliegt. 1: procedure ^^′ ^^ ( ^^ ^^ ) 2 , ^^ ^^ 2: Berechne ^^ = ( ^ ^^ ^^ ^2 )2 3: if ^^ > 1 then 4: Setze κ′ = | ^^| 5: while ^^ > κ′2do 6: Setze κ = κ′ 7: Berechne  κ′ = κ +1 8: Berechne ^^′ ^^zu 9: else 10: Setze ^^′ ^^ = ^^ ^^ Nachfolgend wird die Approximation der transzendenten Funktion ^^(Ϛ) aus der Gleichung (1.3) erläutert. Die Auswertung der in Gleichung (1.3) definierten transzendenten Funktion ^^(Ϛ) ist zu jedem Zeitpunkt tk, zu dem ein neues relatives Winkelinkrement ^^ vorliegt, für ^^ die Berechnung des kumulierten Winkelinkrements ^^ ^^ erforderlich. Im k-ten Zeitschritt ist durch die halbe Norm des Ergebnisses für das kumulierte Winkelinkrement aus dem vorhergehenden Zeitschritt gegeben, also durch = ^^ ^^―1 2 . Die Anwendung des oben zur Renormierung der Winkelinkremente beschriebenen Algorithmus garantiert dabei, dass die Norm ^^ ^^―1 die Bedingung 0 ≤ ^^ ≤ ^^ erfüllt. In dem folgenden Abschnitt werden numerische Approximationen für ^^(Ϛ) hergeleitet, die es erlauben, ^^(Ϛ) effizient und mit der gewünschten Präzision über den gesamten erforderlichen Wertebereich zu berechnen. Die in der Gleichung (1.3) definierte Funktion besitzt die folgende Taylorreihen- entwicklung: ausgewertet ∈ ℕ ausdrücken. Die Riemannsche ist definiert über die Reihe (3.21) Einsetzen der Gleichung (3.20) in die Gleichung (3.19) ergibt dann für die Taylorreihe von ^^(Ϛ) wobei ^^ ^^0+1 das Restglied der Approximation von ^^ durch das Taylorpolynom vom Grad 2 ^^0 in Ϛ ist. Es ist wichtig, dass Ϛ(2n), d. h. die Riemannsche nicht mit dem Argument Ϛ verwechselt wird. Die erfindungsgemäß vorgeschlagene Padé-Approximation ist bei annähernd gleichem Rechenaufwand für die Berechnung erheblich genauer, als die Verwendung des Taylorpolynoms zur Approximation. Die Eignung der Padé- Approximation beruht darauf, dass aus der Definition in Gleichung (3.21) der Riemannschen Mit den ersten denen der Taylorentwicklung übereinstimmen und unter Berücksichtigung der Tatsache, dass nach der Gleichung (3.23) Quotienten aus zwei mit großem Argument in niedrigster Ordnung sich asymptotisch 1 (Eins) annähern, kann der Ausdruck für das Restglied zu allen Ordnungen im Argument Ϛ approximativ als geometrische Reihe resummiert werden. Damit wird das Restglied R erhalten Wenn nun in die Gleichung (3.22) das Restglied durch seine in der Gleichung (3.24) definierte Approximation ^^ ^^0+1(Ϛ) ersetzt wird, dann wird folgende Approximation ^^ ^^0(Ϛ) für ^^(Ϛ) erhalten: Die obige Approximation besitzt nur rationale Koeffizienten, obwohl sie scheinbar transzendente Zahlen wie π und Ϛ(2 ^^) enthält, denn es gilt Ϛ(2 ^^) = ^^( ^^) ^^ 2 ^^, wobei ^^( ^^) ∈ ℚ eine rationale Zahl ist, die von n abhängt. Die ersten in der Gleichung (3.25) auftretenden Koeffizienten lauten explizit (3.26) In Figur 4 sind die relativen Fehler bei Verwendung der Approximation durch die Gleichung (3.25) für verschiedene Werte von ^^0 und bei Verwendung des Taylorpolynoms mit verschiedenem Grad n in Ϛ 2 über den gesamten relevanten Wertebereich 0 ≤ Ϛ ≤ ^ 2^ dargestellt. Es ist wichtig anzumerken, dass für ^^ = ^^0 +2 die beiden Approximationen einen vergleichbaren Rechenaufwand besitzen. Für alle Kombinationen von ^^ = ^^0 +2 ist die Approximation durch die Gleichung (3.25) genauer als die Approximation durch das entsprechende Taylorpolynom. So erreicht die Approximation durch die (3.25) bereits bei ^^0 = 10 über fast den gesamten Wertebereich die in der Software Mathematica voreingestellte Präzision „$MachinePrecision“ = 15.9546 (Zahl der Dezimalstellen für die Darstellung von Machine-Precision Zahlen). Dieses ist mit der Approximation durch Taylorpolynome erst für ^^ = 24 und somit mit ungefähr doppeltem Rechenaufwand möglich. Mit Hilfe des in der Gleichung (3.23) angegebenen asymptotischen Verhaltens der Riemannschen kann der Fehler der Approximation durch Gleichung (3.25) und begründet werden, warum diese eine gegenüber der Approximation durch das entsprechende Taylorpolynom deutlich genauere und damit effizientere Berechnung der Funktion ^^(Ϛ) erlaubt. Die in der Gleichung (3.24) durchgeführte Approximation, auf der die Gleichung (3.25) basiert, verursacht bei Berücksichtigung auch des nächstführenden Terms der asymptotischen Entwicklung (3.23) der Riemannschen folgenden Näherungsfehler Dieser Fehler ist als Funktion von monoton wachsend und damit maximal an der oberen Intervallgrenze, also bei Dort ergibt sich also für den relativen Fehler der Approximation der transzendenten Funktion ^^(Ϛ) durch die Gleichung (3.25) (3.28) Dieser Fehler tendiert für wachsendes ^^0mit exponentieller Geschwindigkeit gegen Null. Dagegen wir der maximale Fehler der Approximation durch das Taylorpolynom vom Grad ^^ = ^^0 +2 in Ϛ 2 durch das Restglied aus der Gleichung (3.22) nach Ersetzung ^^0→ ^^ = ^^0 +2 bestimmt. Für dieses lässt sich für große n nach Ersetzung der Riemannschen durch ihre konstante Asymptotik aus der Gleichung (3.23) eine untere wie folgt angeben Auch wird das Maximum an der Obergrenze des Intervalls von angenommen, also = ^ 2^ bei . Der relative Fehler der Approximation von ^^(Ϛ) durch das Taylorpolynom vom Grad ^^ beträgt daher n zwar gegen es für den Grad ^^ = ^^0 +2 mit vergleichbarem Rechenaufwand ein Fehler von Dieser ist . man einen vergleichbaren Fehler mit der Taylorapproximation erreichen, so findet man durch Gleichsetzen der relativen Fehler die Relation woraus folgt, das ^^ > 2 ^^0 wie der der Approximation (3.25) mittels einer Taylorreihe muss also die Anzahl der zu berechnenden Terme und damit der Rechenaufwand etwa doppelt so groß sein. Die Approximation der transzendenten Funktion ^^(Ϛ) wird nun weiter verbessert. Dabei wird ausgenutzt, dass die Approximation durch die Gleichung (3.25) ein Spezialfall der Padé-Approximation ist, die eine Funktion durch eine rationale Funktion von zwei Polynomen im Funktionsargument approximiert, deren Funktionswert und ersten Ableitungen am Entwicklungspunkt mit denen der zu approximierenden Funktion übereinstimmt. Im Fall der Approximation durch die Gleichung (3.25) ist die Funktion und das Zählerpolynom und das Nennerpolynom besitzen die Grade ^^0 +1 bzw.1 in Ϛ 2. Im Allgemeinen ist die Padé-Approximation von ^^(Ϛ) mit einem Zählerpolynom vom Grad ^^ und einem Nennerpolynom vom Grad ^^ in Ϛ 2 durch gegeben und stimmt am bn Funktionswert B(0) und den ersten ^^ = 1,..., ^^ + ^^ Ableitungen ^^ ^^ ^^ überein. Mit der Zahl ^^0 = ^^0 +2 = ^^ + ^^ zeigt sich, dass bei gegebenen 0 ≤ ^^0 ≤ 10 die Padé-Approximationen mit ^^ = ^^0 ― ^^ und ^^ = ^^ und verschiedenen ^^ = 1,…, ^^0 verschieden genau die Funktion ^^(Ϛ) approximieren. Für die Wahl ^^ = ^^0 ― = ^^ ergibt sich im Intervall 0 ≤ Ϛ ≤ ^^ die Padé-Approximation mit der 2 insgesamt größten Genauigkeit. Diese Spezialfälle bei jedem ^^0 können daher mit der folgenden Gleichung (3.34) für die Padé-Approximation bezeichnet werden Die Koeffizienten lauten wie folgt: Die oben angegebenen ganzzahligen Zähler und Nenner werden durch Linear- kombinationen von Produkten von Ϛ-Funktionswerten erzeugt, wobei jeder Term in der Summe gleiche gerade Transzendentalität besitzt. Von jeder dieser Kombinationen fester Transzendentalität wird dann diese durch Division durch eine entsprechende Potenz von π entfernt. Die sich ergebenden ganzzahligen Zähler und Nenner sollten daher auch in den Beziehungen zwischen auftreten, wie sie z. B. aus den Shuffle- und Stuffle-Relationen folgen. In Figur 5 sind dargestellten relativen Fehler | ^^ ^ bei Verwendung der Approximation von ^^ durch die Gleichung (3.25) für verschiedene Werte von ^^0 (Diagramm a)) und bei Verwendung der Padé-Approximation von ^^ ^^ ^^ ^^é, ^^0(Ϛ) mit Gleichung (3.34) für verschiedene Werte von ^^0 (Diagramm b)) über den gesamten relevanten Wertebereich 0 ≤ ^^ dargestellt. Für ^^0 = ^^ +2 besitzen 2 0 die beiden Approximationen die gleiche Anzahl an Termen. Der Rechenaufwand zur Berechnung von ^^ ^^ ^^ ^^ ^^, ^^0+2(Ϛ) ist sogar im Vergleich zu dem für ^^ ^^ ^^ ^^, ^^0(Ϛ) noch etwas geringer, da wegen der gleichmäßigeren Verteilung des Gesamtgrades ^^0 auf Zähler- und Nennerpolynom weniger Potenzen von Ϛ 2 berechnet werden müssen. Für alle Kombinationen von ^^0 = ^^0 +2 ist die Approximation der Funktion ^^(Ϛ) im Vergleich zu der Approximation durch die Gleichung (3.25) sogar noch genauer und dabei auch noch mit etwas geringerem Rechenaufwand berechenbar als die Approximation durch die Gleichung (3.25). So erreicht die Approximation der Funktion ^^(Ϛ) bereits bei ^^0 = 9 bereits über fast den gesamten Wertebereich die in der Software Mathematica voreingestellte Präzision „$MachinePrecision“ = 15.9546 (Zahl der Dezimalstellen für die Darstellung von Machine-Precision Zahlen), während die Approximation (3.25) erst bei ^^0 = 10 annäherend diese Genauigkeit erreicht. Es sei noch angemerkt, dass anders als der exakte Ausdruck in der Gleichung (1.3) keine der drei Approximationen durch die Gleichungen (3.22), (3.25) und (3.2) einen scheinbaren Pol bei Ϛ = 0 besitzt. Die Approximationen können daher über den gesamten Wertebereich 0 ≤ verwendet werden. Nachfolgend wird die Approximation der transzendenten Funktion ^^ (2) (Ϛ) erläutert. In der Berechnung des kumulierten Winkelinkrementes gemäß Gleichung (1.2) in der in (A.30) gegebenen expliziten Form tritt neben der transzendenten Funktion ^^ (Ϛ) auch die Kombination auf. Diese besitzt wie auch schon ^^ in der Gleichung (1.3) selbst scheinbar einen Pol bei = 0. Allerdings ist die Funktion regulär und bei Ϛ = 0 ergibt sich Analog zum Vorgehen für die im vorhergehenden Abschnitt erläuterte Konstruktion von Approximationen von ^^(Ϛ) können die verschiedenen dort vorgestellten Approximationen für ^^ (2) (Ϛ) konstruiert werden. So ist beispielsweise die Taylorreihe von ^^ (2) (Ϛ) durch gegeben. Wenn die Padé-Approximation von aus der Gleichung (3.34) direkt in der Gleichung (3.36) verwendet wird, liefert dies für diese Approximationen das beste Ergebnis. Dieses ist nicht überraschend, denn die Taylorreihe in der Gleichung (3.37) von ^^ (2) (Ϛ) besitzt eben nicht die für die hohe Präzision der Padé-Approximation verantwortliche Eigenschaft der Taylorreihe von ^^(Ϛ) in der Gleichung (3.22), dass die Koeffizienten für große n exponentiell gegen eine Konstante konvergieren (oder selbst exponentielles Verhalten zeigen). Jedoch ist bei direkter Verwendung von der Gleichung (3.34) in der Approximation nach Gleichung (3.36) der Einfluss des scheinbar vorhandenen Pols noch sichtbar. Von Nachteil ist aber insbesondere, dass bei Verwendung dieser Strategie eine getrennte Behandlung des scheinbar vorhandenen Pols, also des Wertes = 0, notwendig ist. Der zuvor genannte Nachteil lässt sich vermeiden und gleichzeitig lässt sich das beste Approximationsergebnis von allen untersuchten Approximationen konstruieren. Dieses wird erreicht, wenn auf Basis der Padé-Approximation durch Gleichung (3.34) ein Ausdruck für ^^ (2) (Ϛ) berechnet wird, der durch Entfernen des scheinbar vorhandenen Pols regularisiert wird. Hierfür kann die Formulierung der Padé-Approximation (3.34) durch Einführung zweier Polynome wie Die Verwendung dieser Darstellung der Padé-Approximation von ^^(Ϛ) in der Definition von ^^ (2)(Ϛ) in der Gleichung (3.36) erlaubt es, aus der erhaltenen Approximation den scheinbar vorhandenen Pol bei Ϛ = 0 folgendermaßen zu entfernen Hier bezeichnet ^^( ^^ 2 ^^ ) ^^, ^^0(Ϛ) die sich ergebene Approximation, die auf der Padé- Approximation von ^^(Ϛ) basiert aber keine numerisch problematischen scheinbaren Pole mehr enthält und somit regulär ist. In Figur 5 sind die relativen Fehler bei direkter Verwendung der Padé-Approximation (3.25) zur Berechnung von ^^ (2) (Ϛ) und bei Verwendung der regularisierten Approximation ^^( ^^ 2 ^^ ) ^^, ^^0(Ϛ) aus der Gleichung (3.40) für verschiedene Werte von ^^0 über den gesamten relevanten Wertebereich 0 ≤ ^^ dargestellt. 2 Das Diagramm a) zeigt die relativen Fehler wobei durch das Einsetzen der Approximation von ^^ ^^ ^^ ^^ ^^, ^^0 aus der Gleichung (3.34) in die Gleichung (3.36) gegeben ist. Diagramm b) zeigt den relativen Fehler, wenn durch die regularisierte Approximation aus der Gleichung (3.40) gegeben ist. Insbesondere für ^^0 > 8 und kleine ist der regularisierte Ausdruck aus der Gleichung (3.40) erheblich genauer als die Approximation, die auf Einsetzen der Padé-Approximation von ^^(Ϛ) in die Definition der Gleichung (3.36) von ^^ (2) (Ϛ) beruht. Die Ergebnisse der Approximationen können zur Berechnung des kumulierten Winkels aus relativen Winkelinkrementen genutzt werden, um damit die Lage des Objektes zu bestimmen. Hierzu kann die in der Zeile 4 des Algorithmus 1 angegebene Operation verfeinert werden durch 4: Berechne ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ aus ^^ ^ ^^ ^, ^ ^ ^^ ^― ^ 1 und ^^ ^ ^ ^ ^ ^^ ^^ gemäß Gleichung (3.4), d. h. in Gleichung (1.2) mit der expliziten Form nach Gleichungen (3.5a), (3.5). Diese Operation wertet die folgende Gleichung auf Basis der mit Hilfe der erläuterten Padé-Approximationen angepassten Gleichung (3.5) aus: (3.41) Diese beinhaltet nun die folgenden Schritte, wobei angenommen wird, dass ^^ ^^― ^^ entweder direkt vor Aufruf des folgenden Algorithmus 5 in diesem Zeitschritt oder am Ende der Prozedur im letzten Zeitschritt mittels beispielsweise des Algorithmus 2 durch den Repräsentanten mit minimaler Norm ersetzt werden. Algorithmus 5 Berechnung von ^^ ^^ aus ^^ ^^―1 und ∆ ^^ ^^ 1: procedure ^^ ^^ ( ^^ ^^―1 , ∆ ^^ ^^ ) 2: Berechne (θ ^^―1 ) 2 = ^^ ^^―1 ∙ ^^ ^^―1 3: Berechne die Polynome P und aus ^^― ^^ mittels der Gleichung (3.38) mit vorgegebenen N0 4: Berechne ^^ ^^ ^^ ^^ ^^, ^^0(Ϛ) aus P, Q, mittels Gleichung (3.39) 5: Berechne Gleichung (3.39) 6: Berechne ^^ aus ^^ , ∆ ^^ gemäß Gleichung 1.2 unter Verwendung von ^^ ^^ ^^ ^^ ^^, ^^0, ^^( ^^ 2 ^^ ) ^^, ^^0 in Gleichung (3.41) Es wird darauf hingewiesen, dass bei der obigen abkürzenden Schreibweise ^^ ^^ = ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ und ∆ ^^ ^^ = ^^ ^^ ^ ^ ^ ^ ― ^^ 1 ^^ + ^^ zu identifizieren sind. Die Genauigkeit wird durch den Abbruch der Entwicklung der Gleichung (3.41) nach der zweiten Ordnung in ^^ 2 ^^ ^^ ^^ limitiert, was den in der Praxis vor allem relevanten Fall ∆ ^^ ≪ 1 jedoch abdeckt. Bei Bedarf kann eine Erweiterung der Gleichung (3.41) mit den oben beschriebenen Methoden ohne konzeptionelle Hindernisse vorgenommen werden. Im Folgenden wird die Berechnung trigonometrischer Funktionen und der Rotationsmatrix erläutert. Die im Sensorpaket der IMU vorhandenen Beschleunigungssensoren messen je nach Messprinzip entweder direkt die Beschleunigung oder sie messen relative Geschwindigkeitsinkremente, die in Analogie zur Relation (1.1) für die relativen Winkelinkremente durch gegeben sind. Hier bezeichnet die Beschleunigung des Beschleunigungs- sensors bezüglich eines Inertialsystems I, ausgedrückt im Bezugssystem S des Sensors am selben Zeitpunkt t, an dem der Wert der Beschleunigung a(t) vorliegt. Wenn der Sensor also Geschwindigkeitsinkremente misst, so sind diese also aus einer Integration der Beschleunigung im mitbewegten Bezugssystem S entstanden. Insbesondere sind in der Gleichung (1.1) daher auch Drehbewegungen enthalten, die der Sensor während der Integrationszeit ausführt. Die Geschwindigkeits- inkremente gemäß Gleichung (1.1) können daher nicht direkt zur Bestimmung der Geschwindigkeit und Position des Fahrzeugs gegenüber eines Inertialsystems oder anderen Referenz-Bezugssystems verwendet werden. Benötigt werden stattdessen Geschwindigkeitsinkremente, die durch Integration der Beschleunigung in einem zumindest im Integrationsintervall feststehenden Koordinatensystem berechnet worden sind. Das entsprechende Integral ist von der Form wobei die Rotationsmatrix ( ^ ) ^^ die im Sensorsystem S zum Zeitpunkt t gemessene Beschleunigung in ein am Startzeitpunkt ti-1 der Integration dem Sensorsystem entsprechendes aber nicht mit diesem mitrotierendes Koordinaten- system transformiert, d. h. es ist ^ ^ ^ ^ ( ( ^ ^ ^ ^ ^ ) ^― , ^^ 1) ^^( ^^) = ^ ^ ^ ^ ( ( ^ ^ ^ ^ ^ ) ^―1) ^^ ^ ^ ^ ^ ( ( ^ ^ ^ ^ ^ ) ^) , ^^ ^^( ^^′) und somit ist ^ ^ ^ ^ ^ , ^― ^^ 1 ∆ ^^ ^^ das auf das Sensorbezugssystem Si-1 = S(ti-1) bezogene und in diesem auch durch Integration erhaltene relative Geschwindigkeitsinkrement. Um nun eine Anzahl n dieser relativen Geschwindigkeitsinkremente aufaddieren zu können, sind die individuellen Geschwindigkeitsinkremente wiederum in ein gemeinsames Bezugssystem, z. B. SI-1 zum Startzeitpunkt ^^ der Summation, zu transformieren. Die Summation lautet (ohne von Gravitation) also Hierin ist ^ ^ ^ ^ ^ ^ ^ ^ ― ― 1 1+ ^^ ^^ die mit dem vom Startzeitpunkt ^^ ^^―1 aus kumulierten relativen Winkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ 0 ^ ^^ assoziierte Rotationsmatrix. Dieses motiviert, warum die mit dem jeweiligen kumulierten Winkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ 0 ^ ^^ assoziierte Rotationsmatrix ^^ ^ wird. Der Algorithmus ist nicht nur für die oben beschriebene Anwendung, der Verarbeitung der Sensordaten einer IMU und der effizienten Berechnung von Position und Lage eines Fahrzeugs relevant. Darüber hinaus sollte insbesondere die hier vorgestellte effiziente Berechnung trigonometrischer Funktionen vielfältige Anwendungsmöglichkeiten in verschiedenen Bereichen haben, z. B. in Computer Vision. Es wird im Folgenden gezeigt, dass das für das Update des kumulierten Winkelinkrementes ^^ ^^ bereits berechnete Ergebnis für den Funktionswert der Funktion ^^(Ϛ) benutzt werden kann, um effizient und insbesondere ohne Approximation irgendwelcher weiterer transzendenter Funktionen die (trigonometrischen) Funktionen ^^ ^^ 2 ^^ Ϛ2Ϛ und ^^ ^^ ^^2Ϛ zu erhalten. Mit diesen Ergebnissen lässt sich dann direkt die assoziierte Rotationsmatrix ^ ^^^ ^^ ^^ ( ^^) der aktiven Rotation zwischen den Koordinatensystemen F und G berechnen, die mittels der Drehung aufeinander abgebildet werden. Die aktiven Rotation transformiert die Basisvektoren ^^ ^ ^^^ des Koordinatensystems F in die Basisvektoren ^^ ^ ^^^ = ^ ^^^ ^^ ^^(2Ϛ) des nach Anwendung der durch den Winkel ^^ beschriebenen Drehung erhaltenen Koordinatensystems G. Die Kenntnis dieser Matrix ist ein wesentlicher Baustein für die Berechnung der translatorischen Geschwindigkeit und letztendlich der Position aus den von den Beschleunigungssensoren gemessenen Beschleunigungen oder relativen Geschwindigkeitsinkrementen. Es wird Ϛ = identifiziert, wobei θ ist und das erste Argument in Update des kumulierten Winkelinkrementes zugrundeliegenden Gleichung (3.41) bezeichnet und θ ^^―1 = ^^ ^^―1 ∙ ^^ ^^―1 ist und der Repräsentant minimaler Norm des kumulierten aus dem letzten Zeitschritt ist. Allerdings ist der folgende Algorithmus nicht auf die konkrete Anwendung auf die kumulierten Winkelinkremente ^^ ^^―1 beschränkt. Vielmehr könnte ^^1 irgendein Winkel und Ϛ = θ 21 sein halber Betrag sein. Zunächst werden Ausdrücke für (trigonometrische) Funktionen hergeleitet, die von der halben Norm des kumulierten Winkelinkrementes abhängen, also von wobei die Abhängigkeit mittels ^^(Ϛ) ausgedrückt wird. Aus der Definition von ^^(Ϛ) in der Gleichung (1.3) folgt womit nach Einsetzen in die für Argumente 0 ≤ ≤ gültig sind, die folgenden Ausdrücke hergeleitet werden 2 können Das Argument Ϛ = θ 21 = θ ^^―1 ist lediglich die halbe Norm des kumulierten 2 Winkelinkrementes. jedoch die trigonometrische Funktionen der Norm selbst benötigt. Diese können aus den bekannten Relationen zwischen trigonometrischen Funktionen mit einfachem und doppeltem Winkel erhalten werden. Dies führt zu Ein Vorteil der Tatsache, dass die Funktion ^^(Ϛ) nur von der halben Norm des Drehwinkels abhängt, ist also, dass durch die dann notwendige Verdoppelung des Argumentes der trigonometrischen Funktionen die Quadratwurzeln eliminiert werden und sich die Ausdrücke in Gleichung (3.45) ergeben, die rationale Funktionen in Ϛ 2 und ^^(Ϛ) sind. Diese enthalten in Form der Kombinationen ^^ ^^ 2 ^^ Ϛ2Ϛ und ^^ ^^ ^^2Ϛ, die gerade Potenzreihen in Ϛ sind, dann nicht einmal mehr eine Abhängigkeit von der Quadratwurzel aus der Berechnung der Norm selbst. Durch Verwendung der Approximationen gemäß Gleichung (3.25) oder (3.34) können also die trigonometrischen Funktionen durch rein rationale Ausdrücke sehr effizient approximiert werden. Bei Verwendung der Gleichung (3.34) wird In der Figur 6 sind die relativen Fehler der Approximationen von sin ^^ (Diagramm a)) und cos θ (Diagramm b)) mit der Padé-Approximation nach Gleichung (3.34) in den Approximationen gemäß Gleichung (3.46) in dem relevanten Wertebereich 0 ≤ ^^ ≤ ^^ dargestellt. In der Figur 7 sind die relativen Fehler der Approximationen von sin ^^ (Diagramm a)) und cos θ (Diagramm b)) bei der direkten Approximation durch Taylorpolynome vergleichbaren Grades dargestellt. Figur 8 zeigt die relativen Fehler der Approximationen von 2 sin ^ 2^ cos ^ 2^ (Diagramm a)) und 1 ― 2 sin 2 ^ 2^ (Diagramm b)) bei der direkten Approximation durch Taylorpolynome vergleichbaren Grades dargestellt. Es ist zu erkennen, dass der relative Fehler der Approximationen von sin ^^ wegen der Nullstelle bei ^^ = ^^ ansteigt. Die Verwendung der Padé-Approximation durch Gleichung (3.34) in den in diesem Abschnitt hergeleiteten Approximationen gemäß Gleichung (3.46) liefert gegenüber einer direkten Approximation durch ein Taylorpolynom vergleichbaren Grades über den gesamten Wertebereich die genauere Approximation. Die Approximation der trigonometrischen Funktionen mit halbem Winkel durch Taylorpolynome und die sich daran anschließende Berechnung der Funktionen mit dem ganzen Winkel als Argument ist eine deutliche Verbesserung gegenüber der direkten Approximation durch die entsprechenden Taylorpolynome. Allerdings ist die hier entwickelte Methode bei vergleichbarem Rechenaufwand insbesondere für kleinere Winkel erheblich genauer. Ein weiterer Vorteil der Approximation durch die Gleichung (3.46) ist, dass die Bedingung sin 2 ^^ + cos 2 ^^ = 1 unabhängig von der gewählten Ordnung N0 erfüllt wird. Die mit einem vektorwertigen Drehwinkel ^^ assoziierte Rotationsmatrix der aktiven Rotation lautet wobei L ein aus der Matrixdarstellung der drei Lie-Algebra-Generatoren der Drehgruppe SO(3) konstruierter Vektor L = (L1, L2, L3)T ist. Also wird für die Konstruktion der Rotationsmatrix gemäß Gleichung (3.47) insbesondere folgende Kombination benötigt: Damit ergibt sich dann für die Rotationsmatrix gemäß Gleichung (3.47) der folgende Ausdruck Auch dieser Ausdruck ist eine rationale Funktion in und ^^(Ϛ), ∙ ^^ und seine Auswertung erfordert keine Auswertung transzendenter und er enthält nur gerade Potenzen von so dass auch die Quadratwurzel ∙ ^^ nicht berechnet werden muss. Bei Anwendung der Rotation, um den Vektor v in einem anderen Koordinatensystem darzustellen, lässt sich die Gleichung (3.50) wie folgt modifizieren Der Algorithmus 6 berechnet die Rotationsmatrix selbst bei gegebenen ^^, ^^2 und ^^ (Ϛ), z. B. approximiert via Gleichung (3.34). Es ist hier anzumerken, dass darin die Berechnung der Kombination ^^ = nicht erforderlich ist, wenn der Algorithmus 5 zuvor ausgeführt benötigt ebenfalls die Kombination C, so dass das Ergebnis einfach im Rahmen von Algorithmus 5 abgespeichert und dann verwendet werden kann. Algorithmus 6 Berechnung von ^ ^ ^ ^^ ^ ^^ aus ^^, ^^ 2 und B ^ 2^ ^ 1: procedure ^ ^ ^ ^ ^^ ^ ^^, ^^ 2 , B ^ 2^ 2: Berechne ^^ 3: Berechne ^ ^ 2^ aus2 ^ ∙ ^^ 4: Berechne Berechne die Rotationsmatrix aus ^ 2^ ∙ ^^, 2 2 5: ^ 2^ ∙ ^^ , ^^ ^ 2^ mittels der Gleichung (3.50) Der folgende Algorithmus 7 ist eine leichte Modifikation von Algorithmus 6, die zeigt, wie die Koordinaten ^^ ^^ eines Vektors ^^ basierend auf der Gleichung (3.51) direkt im Bezugssystem G aus den Koordinaten ^^ ^^ im Bezugssystem F berechnet werden können, anstatt die Rotationsmatrix zu berechnen. Algorithmus 7 Berechnung von ^^ ^^ aus ^^, ^^ 2 und B ^ 2^ , ^^ ^^ 4: Berechne rechne aus ^^ ^^ ^^ ^^ ∥ ^ 2 5: Be ^^ ^^ ^^ , ^^ , ^^ 2 ^ , ^^ ^^ mittels der Gleichung (3.51) Beide 7, z. B für die effiziente und präzise Berechnung trigonometrischer Funktionen und Drehungen in Echtzeitanwendungen eingesetzt werden, z. B. zum Berechnen von Drehungen im Bereich der Computergrafik / Computer Vision und CAD. Figur 9 zeigt ein Flussdiagramm des Verfahrens, bei dem mindestens ein Gyroskop in x-Richtung, in y-Richtung und in z-Richtung eine Messung der Winkelinkremente ∆ ^^ ^ 1 ^ ,2,3 vornimmt. Diese können aus dem Integral der Winkelgeschwindigkeit ω1,2,3(t) über eine Zeit von einem vorhergehenden Messzeitpunkt tk-1 bis zu einem aktuellen Messzeitpunkt tk gebildet werden. Der Einfachheit halber wurden nur die Gleichungen für die Winkelinkremente dargestellt. Für die Einbeziehung auch der Geschwindigkeitsinkremente ist lediglich jedes Messergebnis eines Gyroskops als um das Messergebnis des Geschwindigkeitsinkrements in der entsprechenden Achsrichtung erweitert zu denken und dann sind korrespondierende Schritte unter Verwendung der Algorithmen für das Geschwindigkeits- und/oder Positionsupdate in analog zu den gezeigten Schritten für die Winkelinkremente durchzuführen. Die Beschreibung der Figur 9 im Folgenden ist der Einfachheit ebenfalls auf die Winkelinkrementupdates beschränkt. Die drei Messwerte für die Winkelinkremente ∆ ^^ ^ 1 ^ ,2,3 können vektorisiert werden, um einen Drehwinkelinkrement-Vektor ∆ ^^ ^^ zu erhalten. Damit können dann die einzelnen Drehwinkelinkremente ∆ ^^ ^^ rekursiv aus dem Exponenten ^^ ^^ ^^ ^^ ^^ mit der Verzögerung (Delay) z-1 zur Einbeziehung des vorhergehenden Drehwinkelinkrement-Vektors ∆ ^^ berechnet werden. Die Berechnung der kumulierten Drehwinkelinkremente ^^ erfolgt mit dem Exponenten ∆ ^^ ^^ , ^ ^ ^ ^ ^^ = 2 ^^ ^ ^^ 2^, ^ ^^^ ^^ ∆ ^^ ^―^ 1 , ∆ 2 ^^ ^^ rekursiv mit einer Verzögerung (Delay) z -1 zur Einbeziehung des vorher berechneten kumulierten Drehwinkelinkrementes ∆ ^ aus den voher bestimmten einzelnen Drehwinkelinkrementen ^^ Das letzte so berechnete kumulierte Drehwinkelinkrement bildet dann den kumulierten Drehwinkel ∆ ^^ ^ ^ ^ ^ ^^ ^^ als Messergebnis.

Claims

Deutsches Zentrum für Luft- und Raumfahrt e.V. Anwaltsakte: Königswinterer Straße 522-524 V/DLR-0799-WO 53227 Bonn-Oberkassel 2022/224 Deutschland Datum: 8. Januar 2024 Patentansprüche 1. Verfahren zum Bestimmen der Lage eines Objektes durch Ermittlung der Dreh- winkel des Objektes, wobei relative Drehwinkelinkremente ∆ ^^ ^^ an aufeinander- folgenden Messzeitpunkten ^^ ^^ mit k = 1, …, n gemessen und für die jeweiligen Messzeitpunkte ^^ ^^ zu kumulierten Drehwinkelinkrementen ^^ ^^ zusammen- gefasst werden, gekennzeichnet durch Berechnen der kumulierten Drehwinkelinkremente ^^ ^^ für die Messzeitpunkte ^^ ^^ iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen ∆ ^^ ^^, die an aufeinanderfolgenden Messzeitpunkten ^ mit k = 1, …, n gemessen wurden, ^ ^^ dem jeweils für das vorhergehende bestimmten kumulierten Dreh- winkelinkrement ^^―1 und d 1 ^^ er transzendenten Funktion ^^(Ϛ) = (1 ― ^^ ^^ ^^ Ϛ) mit als Funktion des kumulierten ^^ ^^―1 zum jeweils vorhergehenden Messzeitpunkt ^^ ^^―1, wobei ein Approximieren der transzendenten Funktion ^^(Ϛ) =Ϛ 1 2(1 ― ^^ ^^ ^^ Ϛ) durch Padé-Approximationen erfolgt. 2. Verfahren nach Anspruch 1, gekennzeichnet durch Berechnen der kumulierten Drehwinkelinkremente ^^ ^^ für einen jeweiligen Messzeitpunkt ^^ ^^ als Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie- Rotationsgruppe SO(3), wobei iterativ kumulierte Drehwinkelinkremente ^^ ^^ für die Messzeitpunkte ^^ ^^ mit k = 1,…, n jeweils aus der Funktion des approximierten Exponenten der Baker-Campbell-Hausdorff-Formel mit jeweils einem zum Messzeitpunkt ^^ ^^ gemessenen relativen Drehwinkelinkrement ∆ ^^ ^^ und dem für den vorhergehenden Zeitpunkt ^^ ^^―1 bestimmten kumulierten JG/JG - 2- Drehwinkelinkrement ^^ ^^―1 oder einem vorgegebenen Start-Drehwinkel- inkrement ^^ 0 als Eingangsgrößen bestimmt werden. 3. Verfahren nach Anspruch 2, gekennzeichnet durch Bestimmen eines kumulierten Drehwinkels ^^ ^^ ^ ^ ^ ^ , ^ ^ ^ ^ ^^ für den Messzeitpunkt ^^ ^^ durch das das zum letzten Zeitpunkt ^^ ^^ des Zeitbereichs I, der die Messzeitpunkte k = 1, …, n umfasst, bestimmte letzte kumulierte Drehwinkelinkrement ^^ ^^ bestimmt ist. 4. Verfahren nach einem der Ansprüche 1 bis 3, dadurch gekennzeichnet, dass das kumulierte Drehwinkelinkrement ^^ ^^ für den Messzeitpunkt ^^ ^ mit der ^ Funktion des approximierten Exponenten ^ ^^ der Baker-Hausdorff- Formel iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen ∆ ^^ ^^, die an aufeinanderfolgenden Messzeitpunkten ^^ ^^ mit k = 1, …, n gemessen wurden, und dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement ^^ ^^―1 berechnet wird, wobei der Beitrag zum Exponenten in nullter Ordnung ^ ^^ ^^ ^ der Beitrag zum Exponenten in erster Ordnung ^^ ^^ ^ ^^ ^^ +  B θ ^^―1 2 ^^ ^^― ^^ 2 ∙ ^^ 2 ^^ ^^ ^^ ^^― ^^ 2 ― θ ^^―1 2 ^^ 2 ^^ ^^ und der Beitrag zum Exponenten in 2   2 ∙ 2 2 × 2 lauten, mit der transzendenten Funktion ^^(Ϛ) = Ϛ2 wobei ein [ ] und die Parameter am und bm von dem Wert ^^0 abhängig vorgegebene Koeffizienten der Padé- Approximation sind, und sich das für den Zeitpunkt ^^ ^^ kumulierte - 3- Drehwinkelinkrement ^^ ^^ jeweils aus der Summe der Beiträge ^ ^^ ^^ 2 ^^ zum Exponenten in der nullten, ersten und mindestens zweiten Ordnung mα = 0, 1 und 2 ergibt. 5. Verfahren nach Anspruch 1 oder 2, dadurch gekennzeichnet, dass die Funktion ^ ^^ die Summe der Argumente der Approximationen in der nullten, ersten, zweiten und dritten Ordnung mα = 0, 1, 2 und 3 umfasst. 6. Verfahren nach einem der vorhergehenden Ansprüche, dadurch gekennzeichnet, dass zur Bestimmung des kumulierten Drehwinkels ΔθI out für den Zeitbereich ^^ ^ ^^ ^^ ^^ den gemessenen Drehwinkel- inkrementen ∆ ^^ bis ^^ für die jeweiligen Zeitpunkte ^^ ^^ mit k = 1,…, n als Eingangsgrößen das Start-Drehwinkelinkrement ^^ ^^ ^ ^ ^ ^ , ^ 0 ^ ^^ für k = 0 mit dem Wert Null vorgegeben wird und die kumulierten Drehwinkelinkremente ^^ ^^ iterativ mit k = 1 bis n aus dem für den vorhergehenden Iterationsschritt k-1 bestimmten oder vorgegebenen kumulierten Drehwinkelinkrement ^^ ^^― ^^ und dem zum Zeitpunkt ^^ ^^, das zum jeweiligen Iterationsschritt k gehörig ist, gemessenen Drehwinkelinkrement ∆ ^^ ^^ als die transzendenten Funktion ^^(Ϛ) = 1 (1 ― ^^ ^^ ^^ Ϛ) umfassenden Funktion des Exponenten der Baker-Campbell- Formel für die Lie-Rotationsgruppe SO(3) mittels einer Padé- Approximation der tranzendenten Funktion ^^(Ϛ) berechnet wird. 7. Verfahren nach einem der vorhergehenden Ansprüche, gekennzeichnet durch Berechnung einer Rotationsmatrix der Rotation zwischen einem ersten Koordinatensystem F des Objektes und einem nach Anwendung einer durch den kumulierten Drehwinkel ^^ ^^ beschriebenen Drehung erhaltenen zweiten Koordinatensystem G mit den Ausdrücken der mittels der Funktion ^^(Ϛ) und insbesondere deren Padé-Approximationen berechneten trigonometrischen Funktionen und Berechnung der translatorischen Geschwindigkeit und/oder - 4- Position des Objektes aus den in das zweite Koordinatensystem transformierten kumulierten Drehwinkeln ^^ ^^. 8. Verfahren zur automatischen rechnergestützten Berechnung von trigonometrischen Funktionen von Winkeln ^^ ^^ oder Drehungen aus gegebenen Winkeln ^^ ^^, gekennzeichnet durch Berechnen der transzendenten Funktion ^^(Ϛ) = Ϛ 1 2(1 ― ^^ ^^ ^^ Ϛ) mit als Funktion des gegebenen Winkels ^^ ^^, wobei ein Approximieren der transzendenten Funktion ^^(Ϛ) = Ϛ 1 2(1 ― ^^ ^^ ^^ Ϛ) durch Padé- Approximationen erfolgt, und Erhalten der trigonometrischen Funktion oder Drehung als Ergebnis der Berechnung der transzendenten Funktion. 9. Verfahren nach Anspruch 8, dadurch gekennzeichnet, dass die Approximation der transzendenten Funktion durch Auswerten der Gleichung ^^ 1 + ∑ ^20 ^^= ^^ ^^Ϛ2 ^^ erfolgt, wobei N0 ein vorgegebener ganzzahliger Wert und die a bm von dem Wert N0 abhängig vorgegebene der Padé-Approximation sind. 10. Verfahren nach Anspruch 9, gekennzeichnet durch Erhalten der trigonometrischen Funktionen als Ergebnis der Padé-Approximation der transzendenten Funktion durch 11. Messeinrichtung zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehraten des Objektes, wobei die Messeinrichtung eine Sensoreinheit zum - 5- Messen von relativen Drehwinkelinkrementen ∆ ^^ ^^ für aufeinanderfolgende Zeitpunkte ^ und eine Datenverarbeitungseinheit hat, die eingerichtet ist, die ^ ^^ an aufeinanderfolgenden Zeitpunkten ^^ ^^ gemessenen relativen Drehwinkel- inkremente ∆ ^^ ^^ zu einem kumulierten Drehwinkelinkrement ^^ ^^ für die jeweiligen Messzeitpunkte ^^ ^^ zu kumulierten Drehwinkelinkrementen ^^ ^^ zusammenzufassen, dadurch gekennzeichnet, dass die Datenverarbeitungseinheit zum Berechnen der kumulierten Drehwinkel- inkremente ^^ ^^ für die Messzeitpunkte ^ iterativ aus einer Anzahl n von ^ ^^ relativen Drehwinkelinkrementen ∆ ^^ ^^, die an aufeinanderfolgenden Messzeitpunkten ^^ ^^ mit k = 1, …, n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement ^^ ^^―1 und der transzendenten Funktion ^^ ^^ ^^ Ϛ) mit als Funktion des kumulierten Drehwinkelinkrementes zum jeweils vorhergehenden Messzeitpunkt ^^ ^^―1 und zum Approximieren der transzendenten Funktion ^^(Ϛ) 1 = ― ^^ ^^ ^^ Ϛ) durch Padé-Approximationen eingerichtet ist. 12. Messeinrichtung zum Bestimmen der Lage eines Objektes und/oder zur automatischen rechnergestützten Berechnung von trigonometrischen Funktionen von Winkeln ^^ ^^ oder Drehungen aus gegebenen Winkeln ^^ ^^ mit einer Datenverarbeitungseinheit, dadurch gekennzeichnet, dass die Datenverarbeitungseinheit zur Durchführung der Schritte der Verfahren nach einem der Ansprüche 1 bis 10 eingerichtet ist. 13. Computerprogramm, umfassend Befehle, die bei der Ausführung des Computerprogramms durch eine Datenverarbeitungseinheit diese veranlassen, die Schritte des Verfahrens nach einem der Ansprüche 1 bis 10 auszuführen.
EP24700394.0A 2023-01-12 2024-01-08 Verfahren und messeinrichtung zur bestimmung der lage eines objektes Pending EP4649284A1 (de)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
DE102023100648.7A DE102023100648B3 (de) 2023-01-12 2023-01-12 Verfahren und Messeinrichtung zur Bestimmung der Lage eines Objekts
PCT/EP2024/050294 WO2024149707A1 (de) 2023-01-12 2024-01-08 Verfahren und messeinrichtung zur bestimmung der lage eines objektes

Publications (1)

Publication Number Publication Date
EP4649284A1 true EP4649284A1 (de) 2025-11-19

Family

ID=89620193

Family Applications (1)

Application Number Title Priority Date Filing Date
EP24700394.0A Pending EP4649284A1 (de) 2023-01-12 2024-01-08 Verfahren und messeinrichtung zur bestimmung der lage eines objektes

Country Status (3)

Country Link
EP (1) EP4649284A1 (de)
DE (1) DE102023100648B3 (de)
WO (1) WO2024149707A1 (de)

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2195811A1 (en) 1996-03-18 1997-09-19 Daniel A. Tazartes Coning compensation in strapdown inertial navigation systems
US5943321A (en) 1997-08-20 1999-08-24 General Datacomm Inc. Circuit set-up and caching for multimedia multipoint servers
US6522992B1 (en) * 2000-05-24 2003-02-18 American Gnc Corporation Core inertial measurement unit
US10274509B1 (en) * 2014-04-09 2019-04-30 Inertialwave Inertial motion tracking device
US10317214B2 (en) 2016-10-25 2019-06-11 Massachusetts Institute Of Technology Inertial odometry with retroactive sensor calibration
DE102022121662B3 (de) * 2022-08-26 2023-09-28 Deutsches Zentrum für Luft- und Raumfahrt e.V. Verfahren und Messeinrichtung zur Bestimmung der Lage eines Objektes

Also Published As

Publication number Publication date
DE102023100648B3 (de) 2024-04-25
WO2024149707A1 (de) 2024-07-18

Similar Documents

Publication Publication Date Title
DE69328850T2 (de) Hybrides Vorwärts-Differenzierungsverfahren und System zur Darstellung von Bezier-Splines-Kurven
EP3329216B1 (de) Bestimmung einer anordnungsinformation für ein fahrzeug
DE69015928T2 (de) Inertial-Transformationsmatrixgenerierung.
DE69228292T2 (de) Gerät zur Erzeugung einer Zwangsbedingung in einem Molekulardynamik-Verfahren
DE112018006161T5 (de) System und Verfahren zum Steuern eines Fahrzeugs
DE102022121662B3 (de) Verfahren und Messeinrichtung zur Bestimmung der Lage eines Objektes
DE4205869A1 (de) Einrichtung zur bestimmung der relativen orientierung eines koerpers
EP0161668A2 (de) Navigationsverfahren für Fahrzeuge insbesondere Landfahrzeuge
DE112023002453T5 (de) System und Verfahren zur Messung der Orientierung eines Objekts auf der Grundlage einer IMU-MARG-Orientierungsschätzung
DE112022007463T5 (de) Schaufelkoordinaten-kalibrierungsverfahren und -einrichtung, aktualisierungsverfahren und -vorrichtung und bagger
WO2015039755A1 (de) Verfahren zur ermittlung einer aktuellen position eines kraftfahrzeugs in einem geodätischen koordinatensystem und kraftfahrzeug
DE2645416A1 (de) Verfahren und anordnung zur ermittlung der verteilung der absorption eines koerpers
DE10312154A1 (de) Verfahren und Vorrichtung zum Ausführen einer Objektverfolgung
DE102023100648B3 (de) Verfahren und Messeinrichtung zur Bestimmung der Lage eines Objekts
Hinze et al. An optimal memory‐reduced procedure for calculating adjoints of the instationary Navier‐Stokes equations
DE69514603T2 (de) Verfahren und Vorrichtung zur Minimierung eines Fehlers infolge einer Störbewegung bei der Bestimmung der Geschwindigkeit in einem Trägheitsmesssystem
DE102018104310A1 (de) Verfahren und System zum Nachverfolgen eines Objekts mittels Doppler-Messinformation
DE60214018T2 (de) Messung geometrischer variablen einer in einem bild enthaltenen struktur
DE102023128626A1 (de) Multimodale Zustandsschätzung mit maskierten Sensormessungen
WO2020160794A1 (de) Verfahren, vorrichtung, computerprogramm und computerprogrammprodukt zum bereitstellen eines bahnverlaufs eines objekts für ein fahrzeug
EP0557592A1 (de) Einrichtung zum Kalibrieren einer Messeinrichtung
DE102005004568A1 (de) Verfahren zur Berücksichtigung von Messwerten von kalibrierten Sensoren in einme Kalmanfilter
DE19617162C2 (de) Verfahren zur Anwendung eines 3D-Gridding-Prozesses in einem Computertomographen sowie Computertomograph zur Durchführung des Verfahrens
Papadakis The planar photogravitational Hill problem
DE102022109438B3 (de) Verfahren und Vorrichtung zur inertialen sensorischen Messung sowie Computerprogramm

Legal Events

Date Code Title Description
STAA Information on the status of an ep patent application or granted ep patent

Free format text: STATUS: UNKNOWN

STAA Information on the status of an ep patent application or granted ep patent

Free format text: STATUS: THE INTERNATIONAL PUBLICATION HAS BEEN MADE

PUAI Public reference made under article 153(3) epc to a published international application that has entered the european phase

Free format text: ORIGINAL CODE: 0009012

STAA Information on the status of an ep patent application or granted ep patent

Free format text: STATUS: REQUEST FOR EXAMINATION WAS MADE

17P Request for examination filed

Effective date: 20250812

AK Designated contracting states

Kind code of ref document: A1

Designated state(s): AL AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HR HU IE IS IT LI LT LU LV MC ME MK MT NL NO PL PT RO RS SE SI SK SM TR

DAV Request for validation of the european patent (deleted)
DAX Request for extension of the european patent (deleted)