Ski Jumping Trajectory Reconstruction Using Wearable Sensors via Extended Rauch-Tung-Striebel Smoother with State Constraints.

To satisfy an increasing demand to reconstruct an athlete's motion for performance analysis, this paper proposes a new method for reconstructing the position and velocity in the context of ski jumping trajectories. Therefore, state-of-the-art wearable sensors, including an inertial measurement unit, a magnetometer, and a GPS logger are used. The method employs an extended Rauch-Tung-Striebel smoother with state constraints to estimate state information offline from recorded raw measurements. In comparison to the classic inertial navigation system and GPS integration solution, the proposed method includes additional geometric shape information of the ski jumping hill, which are modeled as soft constraints and embedded into the estimation framework to improve the position and velocity estimation accuracy. Results for both simulated measurement data and real measurement data demonstrate the effectiveness of the proposed method. Moreover, a comparison between jump lengths obtained from the proposed method and video recordings shows the relative root-mean-square error of the reconstructed jump length is below 1.5 m depicting the accuracy of the algorithm.


Introduction
Acquiring an accurate estimation of the position and velocity of the athlete during a ski jump has always been of great interest to athletes, coaches, and sport researchers to analyze the movement for improving the jumping performance. Traditional solutions to meet this demand are video analysis techniques based on camera recorded videos [1]. However, such a system has the disadvantage that it generally only covers limited parts of the jumping area [2]. To record the entire jump, a large number of cameras would be required making the system expensive and hard to calibrate. Furthermore, camera systems merely provide information about the athlete's position, while velocity and acceleration have to be estimated via numeric differentiation, which may be subject to large error.
In recent years, wearable sensors, such as inertial measurement units (IMUs) and global navigation satellite system (GNSS) receivers, have been widely employed for motion analysis in sports including ski jumping [3]. Chardonnens et al. [4] presented an IMU-based system measuring ski jumping dynamics including the position and velocity of the center of mass perpendicular to the table during take-off. Groh et al. [2] proposed a measurement system with IMUs and a light barrier to estimate the ski velocity and jump length. Blumenbach [5] designed and tested a helmet with a high-precision Global Positioning System (GPS) receiver for the positioning of ski jumpers. Fasel et al. [6] presented a solution to estimate the center of mass and its velocity using differential GNSS and IMUs for alpine skiing. In comparison to video analysis techniques, trajectory estimation using wearable IMUs and GNSS sensors has the advantages of being easy to use and maintain and having full coverage over all phases in a ski jump. For coaches and athletes, such a system provides the possibility to give feedback after each jump regarding the jumper's performance, such as jump length and take-off velocity data.
Trajectory reconstruction techniques are a class of methods to achieve an accurate state estimation of a moving object by properly combining information from the kinematic model with measurement data from sensors, including but not limited to IMUs, magnetometers, and GNSS receivers. Trajectory reconstruction techniques are also termed flight path reconstruction in the aerospace field. A systematic overview of such methods for an aircraft can be found in [7]. Göttlicher and Holzapfel [8] reconstructed an aircraft's flight path using low-cost sensors with the extended Rauch-Tung-Striebel (RTS) smoother. Similar techniques can be also applied to ski jumping applications.
From a methodological point of view, a trajectory reconstruction problem can be interpreted as a state estimation problem: One powerful tool to solve this is the classic RTS smoother, which was first proposed for linear systems in [9] using a linear Kalman filter [10] as the forward-time estimator. The extended RTS smoother is a modification of the classical RTS smoother capable of considering nonlinear systems. It employs the extended Kalman filter (EKF) as the forward-time filter to deal with the nonlinearity of the system. In this paper, the extended RTS smoother is adopted as an estimator. Moreover, we can formulate state equality constraints in a state estimation problem. Different methods to solve state estimation problems with state constraints, which are termed as constrained filtering (or smoothing) problems, are discussed in the survey [11].
In this paper, we present an offline method based on the extended RTS smoother to estimate the trajectory, including the position and velocity, of a ski jumper by the inertial navigation system and GPS integration. Additionally, to fully utilize all available information to improve the estimation quality, the geometric shape of the ski jumping hill is modeled as a set of additional soft constraints and included in the smoothing framework.

Data Acquisition
The measurement data were collected during a summer session in June 2018 on the jumping hill Schattenbergschanze (Hill Size = 106 m) in Oberstdorf, Germany. During this summer session, the hill slope was covered with water-soaked plastic. The measurement setup consists of the following sensors: a Qstarz BT Q1000eX GPS data logger and an iPhone 7, which contains an InvenSense ICM-20600 IMU, and an Alps HSCDTD008A magnetometer. Table 1 summarizes the performance characteristics of the sensors provided by the respective manufacturers. The GPS data logger provides position and velocity measurements with a 10 Hz sample rate. The IMU measures triaxial rotational rates as well as specific forces (i.e., non-gravitational accelerations) at 100 Hz frequency. The magnetometer provides triaxial local magnetic field flux density measurements also at 100 Hz. To the best of our knowledge, the IMU and magnetometer sensor chips in the iPhone 7 provide raw measurement data with a similar accuracy level as previous studies, e.g., [2,4,6,8]. The raw data measured by the smartphone was logged via the application Phyphox [12] and was attached to the jumper's right upper arm using a running armband with the orientation shown in Figure 1. The GPS logger was fixed to the top of the helmet by adhesive tapes for better signal reception.

Methods
In this section, we first introduce the extended RTS smoother framework for a general system in Section 3.1. Next, modifications made for the joint state and parameter estimation as well as the state constraints in Section 3.2. Finally, the system model for ski jumping, which is implemented in the extended RTS smoother framework, is presented in Section 3.3.

Extended Rauch-Tung-Striebel Smoother
The extended RTS smoother is an offline estimation method utilizing the entire batch of measurements over a fixed time interval [16]. It combines a forward filtering pass using the EKF with a backward smoothing pass with the RTS smoother.
A state-space representation for a general nonlinear system described by ordinary differential equations (ODEs) can be written asẋ where x ∈ R n x is the state vector, u ∈ R n u is the known input vector, y ∈ R n y is the model outputs, and z ∈ R n y represents the measurements. Symbols w ∈ R n w and v ∈ R n y denote process and measurement noise vectors. This system can be linearized and discretized into the following form [16]: where H k = ∂g(x k )/∂x. The transformation matrices Φ k−1 , Γ u,k−1 , and Γ k−1 are computed as where , and ∆t k−1 = t k − t k−1 . The noise processes {w k } and {v k } are assumed to be zero-mean, uncorrelated, white Gaussians, i.e.,

Extended Kalman Filter
The EKF is a sequential state estimation method for nonlinear systems [11,16]. First, the initial statex 0|0 and covariance P 0|0 for the EKF are given aŝ wherex 0 is an a priori estimate of the initial states x 0 andP 0 is the covariance of the estimation error ofx 0 . Then, for each time point t k with k ∈ {1, 2, . . . , N}, the EKF consist of a prediction and an update step as follows: 1) Prediction step:x wherex k|k−1 represents the a priori state estimate at time t k based on measurements up to time t k−1 , i.e.,x k|k−1 = E(x k |z 1 , · · · , z k−1 ), and P k|k−1 denotes the covariance of the estimation error ofx k|k−1 .
To evaluate Φ k−1 and Γ k−1 , matrices A k−1 and F k−1 in Equation (7) and Equation (9) are linearized The integral in Equation (12) is numerically solved via the classic 4-th order Runge-Kutta method withũ k−1 = t k −τ 2) Update step: x k|k =x k|k−1 wherex k|k denotes the a posteriori state estimate at time t k based on measurements up to time t k , i.e.,x k|k = E(x k |z 1 , · · · , z k ), and P k|k denotes the covariance of the estimation error ofx k|k . The symbol I n x denotes an n x by n x identity matrix, and H k is approximated by H k = ∂g(x k|k−1 )/∂x .

Rauch-Tung-Striebel Smoother
After the forward-time filtering pass, the EKF estimation of the statesx N|N and its error covariance P N|N is used to initialize the corresponding values in the RTS smoother for the final time point t N : where the superscripts "s" indicate the smoothed estimates. Then, the RTS smoother runs backward in time for t k with k ∈ {N − 1, N − 2, . . . , 0} as where the smoothed estimatex s k can be regarded as the state estimate of time point t k utilizing all measurements, i.e.,x s k = E(x k |z 1 , . . . , z N ) [17], and P s k is the covariance of the estimation error ofx s k .

Adaption for Parameter Estimation and Soft Constraints
To account for unknown parameters and known constraints for the estimation problem, the extended RTS smoother for state estimation is modified by the following adaptations.

Joint State and Parameter Estimation
When the model also includes parameters p ∈ R n p to be estimated, these parameters are treated as additional states withṗ = 0 so that they are compatible with the extended RTS smoother framework. Therefore, with the augmented state vector defined as x a = x p ∈ R n x +n p , the augmented system ODEs becomeẋ

Constrained Filtering
Suppose that a system satisfies the nonlinear state equality constraints as In order to include these into the extended RTS smoother framework, the constraints can be treated as additional pseudo measurements in the form of For the case that equality constraints are satisfied exactly, which is termed as hard constraints in [11,18], the pseudo measurements can be treated as perfect with pseudo measurement noises v c,k = 0 . Contrary, soft constraints are only required to be fulfilled approximately, i.e., c k ≈ h(x k ) . In this case, we assume v c,k are zero-mean Gaussian noises with small covariance Q c , i.e., v c,k ∼ N (0, Q c ) .

We further write Equation (24) as
where z a = z k c k , y a = y k h(x k ) , and v a = v k v c,k . By Equations (22), (25), and (26), we formulated the state-space presentation of an augmented system for state and parameter estimation with state constraints in the form of Equations (1)-(3).

System Model
After introducing the estimation framework for a general system in the previous two subsections, the system model for the ski jumping application is presented in this subsection.

Coordinate Frame Definitions
The first necessary frame is the hill reference frame Ox N y N z N . As shown in Figure 2, its origin O is located at the end of the in-run table. The x N axis and z N axis are both within the symmetric plane of the jumping hill and point horizontally forward and vertically down, respectively. The y N axis is perpendicular to x N and z N axes to form a right-handed coordinate system. r N = [x N y N z N ] denotes the coordinate matrix in the hill reference frame of the jumper's (or, to be exact, the IMU sensor's) relative position vector r with respect to the origin O. Another necessary frame is the local north-east-down (NED) frame Ox O y O z O . The origin of the local NED frame coincides with point O, and the three axes x O , y O and z O point to north, east, and down directions respectively. Finally, the sensor's body frame O B x B y B z B is defined as following: the origin O B is located at the center of the cellphone, and the three axes are fixed to its three axes of symmetry as shown in Figure 3. The coordinates in the NED frame of the relative position vector r are denote by

Conversion of GPS Position Measurements to the Hill Reference Frame
The position measurements from GPS receivers are usually provided in geodetic coordinates, such as the World Geodetic System 1984 (WGS84) frame [19], as a time sequence of latitude µ, longitude λ, and height h. For ski jumping applications, it is more useful to calculate relative positions of the ski jumper with respect to the jumping hill.
With the frame definitions stated before, the GPS-measured positions r WGS,GPS = [µ λ h] in the WGS84 frame are first transformed into the local NED frame r O, GPS = [x O y O z O ] by the following equations: where µ 0 , λ 0 , and h 0 are the latitude, longitude, as well as height of the hill reference frame's origin O, and parameters for the WGS84 reference ellipsoid [19] can be found in Table 2. Table 2. Parameters for the WGS84 reference ellipsoid.

Name Symbol Value
The semi-major axis a 6378137.0 m The first eccentricity e 0.0818191908426 The meridian radius of curvature M µ a 1−e 2 (1−e 2 sin 2 µ 0 ) 3 2 The radius of curvature in the prime vertical N µ a √ 1−e 2 sin 2 µ 0 Then, the GPS-measured positions r N, GPS in the hill reference frame are obtained by where M NO is the transformation matrix from local NED to hill reference frame as where ψ N denotes the hill azimuth angle, which is defined as the rotation angle from the x O axis (the north) to the x N axis, as shown in Figure 2.

Constraint Modeling of the Geometric Shape of a Ski Jumping Hill
The geometric shape of the ski jumping hill is used as additional a priori information and included in the trajectory reconstruction problem as a set of constraints for further improvement of the estimation quality. The mathematical model of a ski jumping hill is described in [20], which was published by the International Ski Federation (Fédération Internationale de Ski, FIS) as the guideline to construct ski jumping hills. Additionally, all necessary model parameters can be found in the jumping hill certificate, which is also issued by the FIS. In our case, we have the certificate of jumping hill -No. 215/GER 29 [21] for the Schattenbergschanze.  Figure 4 shows the geometric profile of a FIS-certified ski jumping hill in the longitudinal plane. The curve from the point A to the origin O is called the in-run where the jumper skis down following the track before take-off. For any point P 1 on the in-run curve, we assume its coordinates to be (x P1 , 0, z IR (x P1 )) where z P1 = z IR (x P1 ) denotes the mapping relationship for the in-run in the vertical plane. Since the trajectory of the jumper before take-off must coincide with the in-run curve, we define the constraints for the in-run curve in the longitudinal and lateral plane as When the jumper is in the in-run area before take-off, i.e., t k ≤ t TO , the constraints should fulfill that c V, IR ≈ 0 and c H, IR ≈ 0 . The curve BC in Figure 4 is the landing area (including the landing slope and the out-run). For a point in the landing area P 2 , its position coordinate in the hill reference frame is ( is the mapping relationship of the landing area in vertical plane. The jumper would land somewhere in between and afterward skiing along this surface. Therefore, we can formulate another constraint as follows: After the jumper's landing, i.e., t k ≥ t TD , the constraints should be active such that c V, LA ≈ 0 . The method for detecting the take-off time point t TO and the touch-down time point t TD using raw measurements is introduced in Appendix A.

Measurement Error Models
To model the measurement error of the inertial sensor, constant bias terms ∆ω and ∆a for the gyroscope and accelerometer as well as zero-mean, white, Gaussian noise terms w gyro and w acc are explicitly considered as where ω B, gyro represents the rotational rates measurement from the gyroscope, a B, acc denotes the measurement from the accelerometer, and ω B as well as a B are the true rotational rates and specific forces respectively. The measurement error model of the magnetometer is considered to be where m B, mag is the magnetometer measurement of the local magnetic field strength, ∆S m represents the vector of the scaling factor error and diag(∆S m ) denotes the corresponding square matrix with S m being its diagonal elements, ∆m denotes the constant bias vector, m B is the true local magnetic field strength, and v mag represents the zero-mean, white, Gaussian measurement noise.
The GPS measurements are modeled as where r N is the true position in the hill reference frame, v O is the true velocity in the local NED frame, r N,GPS and v O,GPS are the GPS measurements accordingly, and v pos and v vel are the measurement noise terms.

State, Output, and Constraint Equations
The system states vector x is given by where v B denotes the triaxial velocity in the body-fixed frame, and q BO represents the attitude quaternions of the body-fixed frame with respect to the local NED frame. The input vector u contains the IMU measurements ω B, gyro and a B, acc as The parameter vector p to be estimated includes the biases of the gyroscope ∆ω, the accelerometer ∆a, the magnetometer ∆m, and the scaling factor error of the magnetometer ∆S m : The process noise vector w consists of influences from both gyroscope and accelerometer: With the definitions above, the state equationsẋ = f (x, u, p, w) stated in Equation (22) can be assembled by the equations of inertial navigation mechanism and kinematics as follows: where M BO is the transformation matrix from the local NED frame to the body-fixed frame as and g O = 0 0 g 0 represents the gravitational vector in the local NED frame. The term a B in Equation (45) where r N can be directly obtained from the state vector, the velocity v O is calculated as and m B, meas is the estimated magnetometer measurement. By taking the magnetometer model Equation (37) into account, m B, meas is estimated by The vector of constraints c is defined to be where c q = q BO 2 − 1 is the constraint for attitude quaternions. Here, c q ensures a unit quaternion, i.e., q BO 2 = 1 . The constraints c V, IR , c H, IR , and c V, LA are geometric constraints given by the shape of the jumping hill, which are introduced in Section 3.3.3.

Results and Discussion
In this Section, trajectory reconstruction results from both simulated and real measurement data are presented. The configuration for the extended RTS smoother is first introduced in Section 4.1. The simulated measurements are used to theoretically validate the purposed method in Section 4.2. Results based on real measurements are presented in Section 4.3. The jump length obtained from both the proposed method and video recordings are compared in Section 4.4 to further validate the results based on real measurement data.

Setting
We first introduce necessary configurations for the extended RTS smoother algorithm. For the initial statesx 0 = r N,0v O,0q 0 , the GPS measurement at initial time t 0 is adopted as an estimation for the initial position and velocity asr N,0 = r N, GPS,0 andv O,0 = v O, GPS,0 . The estimated value for initial attitudeq 0 is calculated via the factored quaternion algorithm proposed in [23] using the accelerometer and magnetometer measurements at time t 0 . The initial guess for all measurement error parameters is chosen to be zero asp 0 = ∆ω 0 ∆a 0 ∆m 0 ∆S m,0 = 0 . For the augmented initial statesx a,0 = x 0p 0 , the estimated covariance matrix for the initial states errorP 0 is set as shown in Table 3. In addition, the covariance matrix of noises Q and R is considered constant and the estimated values of Q and R are listed in Tables 4 and 5 respectively. Within the model, the hill azimuth angle ψ N is set as 312.6 • for the jumping hill Schattenbergschanze. The gravitational constant at the site is calculated to be g 0 = 9.8053 m/s 2 according to the Earth Gravitational Model 1996 [19]. Table 3. Diagonal elements of the estimated initial states error covariance matrixP 0 .

System States Symbol Estimated Variance
Position Table 4. Diagonal elements of the estimated process noise covariance matrix Q.

System Inputs Symbol Estimated Variance
Gyroscope

Validation by Simulated Measurement Data
To validate the accuracy of the proposed algorithm, an artificial measurement data set is generated by assuming a true value and adding sensor error and noise. This data set is generated based on the result from real measurements as a reference trajectory, ensuring similar observability of the estimation problem as in the real world. Trajectory optimization techniques with least-squares costs are used to obtain a trajectory with similar acceleration and angular velocities as the measured data. By adding path constraints to the trajectory optimization, the true trajectory also fulfills the constraints in Equations (32)-(34). The true value for the generated trajectory also complies with the kinematic models in Equations (44) Tables 4 and 5 to generate artificial measurements.
The trajectory is reconstructed with only the noisy artificial measurements known to the extended RTS smoother. The reconstruction result is presented in Figure 5. Here, it can be observed that the smoothed trajectory agrees closely with the true reference. To quantify the error of the estimated position and velocity, we calculated the error by where variables with subscript "ref" are the true values from the generated data. The error of these estimated positions and velocities are presented in Figure 6a,b respectively. To intuitively display the attitude, the quaternions q BO are converted to the corresponding Euler angles (in z-y-x order): the bank angle φ, the pitch angle θ, and the azimuth angle ψ by (56) Figure 5. Comparison of the extended RTS smoother reconstructed trajectory (color-coded line), the generated reference trajectory (black line), and the GPS measured trajectory (purple dots) plotted on the ski jumping hill model.
Then, we use the error in Euler angles to represent the attitude estimation error as The Euler angle estimation error is presented in Figure 6c. The root-mean-square (RMS) error for position, velocities, and Euler angles are summarized in Table 6. This result shows that the proposed algorithm can successfully reconstruct the trajectory and that the position accuracy for the assumed noise covariance is in decimeter magnitudes.

Validation by Real Measurement Data
Now, we apply the extended RTS smoother on real measurement data to the ski jumping trajectory reconstruction problem. The results for a reconstructed ski jump are presented in the following part.
In Figure 7, the three-dimensional trajectories of the jumper are plotted with a digital ski jumping hill model based on [20,21]. The purple solid line with dots shows the converted raw GPS position measurements r N, GPS where each dot on the line represents a measurement sample. Furthermore, the line consisting of color-coded dots represents position reconstructed by the extended RTS smoother, where the color of one dot represent the velocity Note that we use the relative local height h N here instead of z N in the figure for convenience. The relationship between h N and z N is directly h N = −z N . The subplot on the top-right corner of Figure 7 presents the projection of the jumping trajectories on the x N -z N plane (focusing mainly on the flight part). The comparison between GPS-measured and extended RTS smoother reconstructed position as well as velocity time histories is presented in Figures 8 and 9 respectively.
Both GPS measurements and the smoother estimation in Figure 7 show the trajectory of an entire jump from start to stop. By comparing the GPS-measured trajectory to the jumping hill model for the in-run and out-run parts, it shows that the GPS height channel contains relatively large errors. After the jumper starts his run along the track during the in-run phase, the GPS-measured trajectory indicates that the jumper starts to move horizontally, but remains at the same height vertically. For part after the touch-down in the landing area, the GPS-measured trajectory shows large error leading to the position being above the ground. Therefore, in this case, the GPS-measured trajectory is not adequate for the demand for analyzing the ski jump limited to its accuracy. It is worth mentioning that the GPS measured trajectory does not necessarily encounter the same type of error in the height channel as in our case. Nevertheless, due to the trilateration positioning principle of the GPS technology employed, the GPS height measurement is far less accurate than the latitude and longitude measurements. The GPS service standard [24] reported that the GPS has a global average positioning accuracy of ≤ 9 m (95%) horizontal error and ≤ 15 m (95%) vertical error.
On the other hand, besides the GPS measurements as a data source, the extended RTS smoother uses IMU and magnetic measurements as well as geometric information of the jumping hill to perform data fusion. The application of the RTS smoother reconstructs the trajectory despite a lower GPS measurement frequency and the existence of some GPS height errors as shown in Figures 7-9. In Figure 7, the reconstructed trajectory shows the entire process of the jump. The jumper first slides down along the in-run track with increasing speed V S . For the flight part, this can be clearly observed from both the three-dimensional plot and the subplot in Figure 7. It can be observed that the height above the ground is decreasing until the jumper's touch-down with the speed further increasing. After the touch-down, the jumper skis down on the landing area of the hill with his speed still accelerating. After reaching the flat part at the hill bottom, the jumper decelerates due to the friction and his break action and then comes to a full stop. In comparison to the GPS measurement, the estimated position and velocity are improved especially in the height channel, which reveals more information for the jump analysis, and the sampling frequency on the trajectory is also increased.
In Figure 8, we can observe that two lines fit well for x N . For y N in the extended RTS smoother result (the green line), the constraint c H, IR correct the position y N to 0 before take-off. Afterward, the reconstructed trajectory matches the GPS measurements with some differences. For the h N channel, as also mentioned before, the raw GPS measurement has a relatively large error which can be clearly observed in Figure 7. The smoother reconstructed height by considering the vertical constraints c V,IR and c V,LA is moving in the in-run curve before the take-off and on the surface of the landing area after the touch down (see Figure 7). The height history is closer to a real physical movement. In Figure 9, for horizontal velocity v O, x and v O, y , the reconstruction and measurement fit accurately. The extended RTS smoother result gives a bit different estimation of v O,z compared to the GPS measured one. Considering the fact that after the jumper starts gliding down around 610.5 s, the vertical velocity afterward should be positive (moving downwards). So, the vertical speed measured by GPS also contains an error for the in-run part which is associated with the height error. Since the extended RTS smoother considers kineticsḣ N = −v O,z , the smoother estimated vertical speed is linked to the height estimation. Combining the result in Figures 7 and 8 for the in-run part (around 610.5 s to 616.8 s), the positive v O,z estimates the smoother giving out is more close to the fact that jumper is moving vertically down in the in-run. This indicates that the extended RTS smoother improves the velocity comparing to the GPS with the assistance of information from other sources.    The improvement of the estimation accuracy is achieved by fusing data from different sources including measurement sensors, kinematics, and geometric information of the hill. The gyroscope and accelerometer data, after removing the estimated bias error, are used for propagating the states including attitude, velocity, and position together with a priori knowledge of the uncertainties. The GPS measurements and the magnetometer are included as additional data sources for the correction step. Furthermore, the introduction of the geometric constraints from the jumping hill model provides an additional information source, which contributes to the estimation. The constraints c V, IR and c V, LA , as shown in Figure 10, establish a strong relationship between states x N and z N on the vertically plane for each part of the trajectory (before take-off and after touchdown).
Although the major aim of applying the extended RTS smoother is to obtain a better position and velocity estimation, it is worth to mention that the RTS smoother also provides an attitude estimation of the IMU sensor. Furthermore, the estimation of the IMU attitude is essential to attribute the accelerometer measured specific forces into the correct direction (see Equations (45) and (49) ). By applying Equations (54)-(56) , the converted Euler angles reconstructed by the extended RTS smoother are shown in Figure 11. From Equation (37), we know that the magnetic field measurements change with the attitude. Therefore, the similar tendency of extended RTS smoother estimated m B, meas and magnetometer measurements m B,mag in Figure 12 indicates good attitude estimation.
To conclude the result presentation, measurements from the gyroscope and the accelerometer are shown in Figures 13 and 14.

Validation by Jump Length
To validate the trajectory reconstruction results, we calculate the jump length results for the five collected trials obtained from the trajectory reconstruction L TR for comparison with the jump length measurements from video recordings L VR . The process of determining the jump length from the video recordings is explained in Appendix B. The validation result of five different jumps is shown in Table 7. Table 7. Comparison between jump length obtained from the trajectory reconstruction results (proposed method) and video recordings (reference). The error ∆L in Table 7 is defined as ∆L = L TR − L VR , and the RMS error for the jump length estimation is therefore ∆L = 1 5 ∑ 5 i=1 ∆L 2 = 1.0 m. Considering that the estimate of the L VR value could contain an error of 0.5 m as worst case, we calculate the upper bound of the RMS error as ∆L ub = 1 5 ∑ 5 i=1 (|∆L| + 0.5) 2 = 1.5 m .

Conclusions and Discussion
This paper presents a novel method for trajectory reconstruction in ski jumping using wearable sensors by applying the extended RTS smoother with state constraints. The result of simulated measurement data validates the proposed method and shows that it has a theoretical accuracy of decimeter magnitude in positions. The results based on real measurements show the method can successfully reconstruct the trajectories of ski jumps. By formulating the geometric information of the jumping hill as soft equality constraints, the estimation accuracy of positions and velocities is largely improved comparing to the GPS measurements. Jump length data for the five collected trails obtained from video recordings demonstrate good accordance with the reconstruction results with the root-mean-square error upper bound of 1.5 m. This provides strong supporting evidence for the validity of the proposed method.
For future research, a quantitative assessment of the accuracy of the proposed method by further real measurement data needs to be carried out. Therefore, another measuring system providing higher accuracy of position and velocity, e.g., a well-calibrated camera system, needs to be employed in order to compare the results of the two systems. For future device developments, wearable sensors providing better position and velocity measurement accuracy, such as differential GPS receivers, which could provide centimeter-level accuracy [6], are recommended to gain better estimation of the jumper's trajectories. In addition, GPS receivers with higher sampling frequency can be tested which could potentially improve the results. On the other hand, due to the importance of size and weight in the performance of ski jumping, the wearability of such sensors is still a major aspect to be considered. Future measurements should include more trials performed by different athletes in order to increase the statistical confidence. After further development, this technique is promising to be implemented to provide quantified motion feedback in ski jumping training. in Figure A1. Then, two peaks (local maximums) with very large a norm and time difference more than 1.5 s (based on experience) are selected by the following rule as shown in Figure A1 as circles: (a norm (t k ) | V GPS (t k ) > 10 m/s, −1.5 s ≤ t k − t peak1 ≤ 1.5 s) .
These two peaks indicate large ground reaction forces for take-off and landing. Finally, the take-off time points t TO is determined as the earlier peak. As presented in [26,27], the actual touch-down point in the landing process is a short moment before the maximal ground reaction peak, approximately 0.05 s as we estimated. Therefore, the touch-down time points t TD is determined as 0.05 s before the later peak as:  Figure A1. An illustrative example on detecting the take-off and touch-down time points by the raw measurement data.

Appendix B. Jump Length Detected by the Video Recordings
In this section, the process for detecting the jump length detected using the video recordings is explained. Cameras with a frame rate of 50 fps are set up at the side of the landing area to record the jumper's landing movement. Figure A2 shows a picture of overlaid snapshots from one of the video recordings. The time difference between each pair of neighboring snapshots is 0.1 s (with a distance of 5 frames). First, the frame for the jumper's touchdown determined indicated by both skis basically touching the ground. Then, the jump length is detected by the middle point of the ski with the help of the distance markers on the other side of the hill, as shown by the red dash lines in Figure A2. Overall, we consider the detected jump length is within an accuracy of 0.5 m. Figure A2. A picture of overlaid snapshots from video recordings to illustrate the jump length detection from video recordings.