Attitude Estimation of Underwater Vehicles Using Field Measurements and Bias Compensation

This paper proposes a method of estimating the attitude of an underwater vehicle. The proposed method uses two field measurements, namely, a gravitational field and a magnetic field represented in terms of vectors in three-dimensional space. In many existing methods that convert the measured field vectors into Euler angles, the yaw accuracy is affected by the uncertainty of the gravitational measurement and by the uncertainty of the magnetic field measurement. Additionally, previous methods have used the magnetic field measurement under the assumption that the magnetic field has only a horizontal component. The proposed method utilizes all field measurement components as they are, without converting them into Euler angles. The bias in the measured magnetic field vector is estimated and compensated to take full advantage of all measured field vector components. Because the proposed method deals with the measured field independently, uncertainties in the measured vectors affect the attitude estimation separately without adding up. The proposed method was tested by conducting navigation experiments with an unmanned underwater vehicle inside test tanks. The results were compared with those obtained by other methods, wherein the Euler angles converted from the measured field vectors were used as measurements.


Introduction
Navigation is one of the fundamental technologies for underwater vehicles and robots. The attitude and velocity of a vehicle is required for the navigation and control of the vehicle [1,2]. This paper proposes a method that improves the accuracy of attitude estimation for which measurements sensed by the Microelectromechanical Systems-Attitude and Heading Reference System (MEMS-AHRS) are used. The MEMS-AHRS measures the acceleration, magnetic field, and angular rate of the sensor in three directions. The measured acceleration is regarded as the measurement of the gravitational field, and the measured magnetic field is regarded as the measurement of the geomagnetic field specific to a geographic location. The MEMS-AHRS has a small volume and is an inexpensive, lightweight, and convenient sensor for determining attitude. In this paper, MEMS-AHRS will be referred to simply as attitude and heading reference system (AHRS).
Amongst the three measurements by an AHRS, the magnetic field measurement is highly susceptible to distortion, owing to the hard-iron and soft-iron effects. Additionally, the scale factor error and non-orthogonality of the sensor construction make the measurement less reliable. The uncertainty In these applications, PF outperforms other methods in dealing with highly nonlinear and complicated measurement model where explicit mathematical description of the environment is impractical.
The complementary filter (CF) method fuses the complementary spectral characteristics of the field measurements and angular rate measurement of the AHRS for attitude estimation [14,15]. Although the angular rate measurement yields an attitude change that is accurate and useful for short-term use, the attitude calculated from the change exhibits drift with time. Thus, its long-term reliability is compromised. On the contrary, the roll and pitch calculated from the acceleration and magnetic field measurement provide the absolute attitude, which is robust in the long term without drift. However, the roll and pitch are not suitable for determining the relative change of attitude in the short-term. Based on this complementary property, the attitude calculated from the angular rate is passed through a high-pass filter, while the attitude calculated from the acceleration and magnetic field is passed through a low-pass filter. Then, the filtered attitudes are combined by the CF, which is being further developed to adjust the gains based on the measurement of acceleration for high acceleration applications [16]. The nonlinear explicit complementary filter (NECF) is partially similar to the CF and has been widely used with many modifications [2,17].
The NECF was derived as an observer in a special orthogonal group SO(3) by exploiting the geometry of the SO(3) [15,17]. It uses the difference between the detected field measurements and the presumably true field through the cross product of the measured field and the true field pair. The method estimates the rotation matrix by reflecting the attitude error on the bias in the angular rate, while the bias in the measured field is not explicitly addressed. The NECF uses only the horizontal component of the magnetic field measurement. In the context of an invariant observer [18], it has been demonstrated that a particular choice of gains reduces the invariant observer to NECF [19].
The problem of representing the attitude and rotation has been investigated extensively because the attitude and rotation constitute SO(3), which is not a linear space. The methods of invariant EKF (IEKF) systematically resolve the nonlinear system estimation problem by using a transformation group based on Lie algebra [19][20][21]. The application of IEKF depends on the measurements available to the specific problem. Therefore, it is still hard to apply IEKF to a general attitude estimation problem, where the AHRS measurements and the additional field measurements are utilized.
The sine rotation vector (SRV) method has been proposed to relieve the problem arising from the Euler angle representation of measurement innovation [22]. As a measurement innovation for the KF application, the SRV method replaces the algebraic difference between the measured Euler angles and the predicted Euler angles with the sine rotation vector. Although the SRV method improves the estimation performance, it still requires the conversion of field measurements into Euler angles, and also requires additional computation to calculate the SRV. The EKF, CF, and SRV basically use the Euler angles converted from the acceleration and magnetic field measurement. Therefore, they all suffer from a distorted magnetic field measurement and uncertainty in the yaw calculation.
This paper describes a method that overcomes the problem of distortion-prone magnetic field measurement and accumulated error in yaw calculation. The proposed method fully utilizes the field measurements and calculates yaw independently from the roll and pitch by using all components of the measured acceleration and magnetic field directly in the measurement update stage without converting them to Euler angles. Additionally, the proposed method estimates and compensates for the magnetic field bias, and uses the magnetic field vector to estimate the roll and pitch as well as the yaw. The proposed method uses the well-established EKF method, which is easily expandable when additional field measurements other than gravity and magnetic field are available. For example, if a landmark is recognized, then the vector directional to the landmark is regarded as the measured field. Thus, a simple augmentation of the measurement variables and measurement model equations is carried out by appending the additional field, and this completes the EKF adaptation to this particular problem.
Section 2 formulates the problem in a manner suitable to the use of field measurements. Section 3 explains the proposed approach in the order of algorithm flow: prediction, bias estimation and compensation, application of measurement model to field measurement, and correction. This section also discusses an optional alternative, which can be used when the vertical component of the magnetic field bias is poorly compensated and the incorrect compensation degrades the estimation performance. Section 4 describes the experiments using unmanned underwater vehicles (UUVs) navigating in test tanks. The results compare the attitude estimation performance and trajectory estimation performance for which the attitude estimated by the proposed method is used. Finally, Section 5 concludes this paper along with suggestions for future work related to the proposed method.

Nomenclature
The notations used for the derivation of the method are as follows: attitude of a robot at time t; , and ψ(t) indicate the roll, pitch, and yaw, respectivelŷ x − (t) attitude predicted for time t, before being corrected using the measurements; attitude estimated for time t through the prediction and correction procedure; T a(t) acceleration measured at time t in the instrument coordinate frame; m(t) magnetic field measured in the instrument coordinate frame; m g (t) geomagnetic field represented in the North-East-Down (NED) coordinate frame appropriate for the geographical location; m g (t) = (m g,x (t), m g,y (t), m g,z (t)) T z u (t) normalized field measurement at time t; z u (t) = (a u (t), m u (t)) T u(t) linear velocity measured by a Doppler velocity log (DVL) in the instrument coordinate frame; u(t) = (u(t), v(t), w(t)) T ω(t) rotational (angular) velocity measured by a gyroscope in the instrument coordinate frame; ω(t) = (p(t), q(t), r(t)) T f(·) motion function relating attitude x(t) and angular velocity ω(t) in the instrument coordinate frame to the robot's angular velocityẋ(t);ẋ(t) = f (x(t), ω(t)) h(·) measurement function relating state x(t) to measurement z u (t); z u (t) = h(x(t)) t k the kth discretized sampling instant of time ∆t k time period between t = t k−1 and In many studies, an underwater robot or vehicle is often referred to as a UUV. Thus, in the following sections, the terms UUV, vehicle, and robot will be used interchangeably.

Problem Formulation
The problem being considered in this paper is described as follows.
Attitude Estimation Problem: Find the attitudex(t) of an underwater robot at time t = t k , given the measurements a(t k ) andm(t k ), and attitudex(t k−1 ) estimated at time t = t k−1 , linear velocity u(t k ), and angular velocity ω(t k ).
The state x(t) to be estimated consists of the Euler angles of a robot; namely, the roll φ(t), pitch θ(t), and yaw ψ(t).
The acceleration measurement a(t) = (a x (t), a y (t), a z (t)) T is normalized to a u (t) = (a u x (t), a u y (t), a u z (t)) T . The magnetic field measurementm(t) = (m x (t),m y (t),m z (t)) T is calibrated to m(t) by subtracting the bias from the measurement. Then, the bias calibrated magnetic field measurement m(t) is normalized to m u (t) = (m u x (t), m u y (t), m u z (t)) T . The normalized measurements a u (t) and m u (t) are merged to form the measurement vector z(t), as expressed by Equation (2): The field vectors a(t) andm(t) are measured in the instrument coordinate frame. It is assumed that the vehicle coordinate frame coincides with the instrument coordinate frame. In the development, the reference coordinate frame, which is otherwise called the world coordinate frame, indicates the North-East-Down (NED) coordinate frame [23].
A vehicle's linear velocity u(t) = (u(t), v(t), w(t)) T is detected by a DVL, and the detected velocity is used to dead-reckon the vehicle location that will be used to verify the applicability and performance of the proposed attitude estimation method. The rotational velocity of the vehicle ω(t) = (p(t), q(t), r(t)) T is detected by a gyroscope to provide the time update of the state. The linear velocity u(t) and rotational velocity ω(t) are also measured in the vehicle coordinate frame, which is adjusted to coincide with the instrument coordinate frame.
The procedure for attitude estimation consists of prediction, bias estimation, bias compensation, and correction, and will be described in the following sections.

Attitude Estimation
The attitude estimation process consists of the prediction and correction of the attitude and attitude covariance. The process uses a magnetic field measurement, whose bias is estimated and compensated. This section explains the prediction stage (Section 3.1), procedure for the estimation and compensation of the magnetic field bias (Section 3.2), derivation of measurement model (Section 3.3), and correction stage (Section 3.4). Additionally, an optional method using only the horizontal component of the magnetic field is applicable when the severe distortion of the vertical component of the magnetic field measurement degrades the attitude estimation. This method is presented in Section 3.5.

Prediction of State and State Covariance
The time update, which is otherwise called the prediction or propagation step, predicts the state for time t k asx − (t k ) and the covariance P − (t k ) of the predicted statex − (t k ). The prediction step updates the previous estimatesx(t k−1 ) and the covariance P(t k−1 ) by using the differential equation for attitude transition. The differential equation model of the attitude is expressed by the following Equation (3): The time derivativeẋ(t) of the state in Equation (4) is a function of the measured angular velocity ω(t) = (p(t), q(t), r(t)) T and state x(t) = (φ(t), θ(t), ψ(t)) T . In these equations, S(·), C(·) and T(·) represent sin(·), cos(·), and tan(·), respectively. In this model, v(t), which is assumed to be an additive and Gaussian with zero mean, is the noisy uncertainty involved in the state transition. The prediction of the state is based on the derivative ofx(t) as expressed by Equation (5): In Equation (6), the derivativeṖ(t) is used to predict the state covariance P(t): where, The Jacobian matrix F(t) is expressed by Equation (8): (8) where the time variable (t) is deleted from p(t), q(t), r(t),φ(t), andθ(t), for notational simplicity.
Using Equations (5) and (6), the state and covariance are predicted for time t k from statex(t k−1 ) and covariance P(t k−1 ), which are estimated at time t k−1 . The predicted statex − (t k ) and covariance P − (t k ) can be calculated by the numerical integration of (5) and (6) for time t ∈ [t k−1 , t k ). One of the simplest methods to predict them is to use Equations (9) and (10): This method can be used under the provision that ∆t k = t k − t k−1 is small enough for derivativeṡx (t k−1 ) andṖ(t k−1 ) to approximate the derivativesẋ(t) andṖ(t), respectively, for the time period t ∈ [t k−1 , t k ).

Bias Estimation and Compensation of Magnetic Field Measurement
After the prediction procedure described in Section 3.1, the predicted attitude and covariance will be corrected by comparing the measurements with the predicted measurements which are calculated based on the predicted state values. However, the magnetic field measurement is corrupted by noise and bias. Therefore, before applying the measurement update procedure, the bias changing with time and the location must be estimated and compensated for. The proposed method estimates the bias by a Kalman filter (KF) using the angular rate measurement [24].
The measurement of the magnetic fieldm(t) includes bias b m (t) and Gaussian noise n m (t), as modeled in Equation (11): where r m tr (t) is the true magnetic field in the reference coordinate frame, whilem(t) represents the measured magnetic field in the instrument coordinate frame. Ignoring the Gaussian noise n m (t) results in the true magnetic field r m tr (t) in the fixed world coordinate frame of Equation (12): In Equation (12), r v R(t) is the rotation matrix that transforms a vector in the vehicle coordinate frame to the vector in the fixed world coordinate frame, expressed by Equation (13): Differentiating Equation (12) results in Equation (14): When deriving Equation (14), it is assumed that the true magnetic field at a location does not change rapidly. Thus, there is no substantial change in the magnetic field with time. Additionally, the same assumption applies to the bias b m (t). Applying the relationship r vṘ (t) = r v R(t)S(ω(t)) to Equation (14) gives the following equation: where the skew-symmetric cross product matrix S(ω(t)) is defined as Equation (16): Multiplying r v R −1 (t) on both sides of Equation (15) yields the differential equation of the magnetic field measurementm(t):m Equation (17) can be rearranged to form the state transition Equation (18) for the KF application to estimate the b m (t) bias along with measurementm(t): The measurementm(t) is obtained from (m(t) b m (t)) T as Equation (19): In the derivation of Equation (18), it is again assumed that bias b m (t) does not change substantially with time, as assumed in the derivation of Equation (15). Equations (18) and (19) correspond to the process and measurement model in the application of the KF to the estimation of the bias b m (t) [24]. Applying the typical KF procedure using the models represented by Equations (18) and (19) provides an estimate of the bias in the magnetic field measurement [23]. The bias calibrated magnetic field expressed by Equation (20) is used in the remainder of the attitude estimation procedure: where b m (t) is provided by the KF, which is implemented using Equations (18) and (19).

Measurement Model and Linearization of the Measurement Model
The proposed method uses the measured magnetic field and acceleration for the correction of the predicted statex − (t k ), and the state covariance P − (t k ). The measurement model is the function which maps the state onto the magnetic field and acceleration measurements. The measurement model is nonlinear and needs to be linearized to apply the EKF approach.
First, the measurement model function h(x(t)) is derived. For the attitude estimation, the normalized measurement z u (t) is used as the measurement z(t): The derivation of the measurement function for the acceleration is based on the assumption that the acceleration of the vehicle motion is negligible in comparison with gravity. The assumption ascertains that the normalized acceleration vector a u (t) = (a u x (t), a u y (t), a u z (t)) T corresponds to the unit vector in the upward direction, which is (0, 0, −1) T in the NED coordinate frame. Therefore, the following Equation (22) describes the relationship between the unit acceleration vector and the attitude: Equation (22) describes the measurement model for the acceleration measurement. The measurement model for the magnetic field measurement is based on the relationship between the geomagnetic field and the measured magnetic field. The measured magnetic field is the geomagnetic field in the sensor coordinate frame. The geomagnetic field vector m g consists of the northerly intensity m g,x , easterly intensity m g,y , and vertical intensity m g,z (positive downwards) of the geomagnetic field at a certain location on Earth: The geomagnetic field value depends on the location and altitude on Earth. In this paper, the latest model of the world magnetic model for 2015-2020 is used [25]. The magnetic field measurement m(t) is determined by applying rotation to the geomagnetic field m g , as expressed by Equation (24): The magnetic field m(t) of Equation (24) is normalized to m u (t) to constitute the measurement vector z, as expressed by the following Equation (25): where m u g = (m u g,x , m u g,y , m u g,z ) T = m g m g = m g m g,x 2 + m g,y 2 + m g,z 2 .
In Equation (25), m g is normalized to m u g . Using Equation (13), m u x (t), m u y (t), and m u z (t) are derived as follows: Using Equations (22) and (25), the measurement model is derived as the function h(x(t)) expressed by Equation (29): The linearized measurement model H(t k ) is derived from Equations (26)- (29): In Equation (30), the time index (t k ) is deleted for notational efficiency. The partial derivative elements of matrix H(t k ) are derived as follows: The linearized matrix H(t k ) is used in the correction procedure described in Section 3.4.

Correction of Predicted State and State Covariance
The state representing the attitude is finally estimated asx(t k ) by adjusting the predicted statê x − (t k ) as expressed by the following Equation (43): The Kalman gain K(t k ) in Equation (43) is determined from the linearized measurement model H(t k ) and the predicted covariance matrix P − (t k ): The predicted covariance P − (t k ) is also corrected to P(t k ) by Equation (45): Equations (43) and (45) complete the attitude estimation procedure. The estimated statex(t k ) and covariance P(t k ) of Equations (43) and (45) are used in the next iteration of the EKF procedure to predictx − (t k+1 ) and P − (t k+1 ) using Equations (9) and (10).

Use of Horizontal Component of Magnetic Field Measurement
In the proposed field measurement approach, the magnetic field measurement affects the roll and pitch estimation, and the yaw estimation. The acceleration measurement plays a major role in estimating the roll and pitch, and the estimated roll and pitch have enough accuracy and stability for use in underwater navigation, where the acceleration of vehicle motion is negligible in comparison with gravity. Because the magnetic field measurement is more vulnerable to environmental distortion than the acceleration measurement, the use of magnetic field measurement occasionally deteriorates the estimation of roll and pitch if the magnetic field measurement is not effectively calibrated to reduce the distortion effect. This happens when the range of the roll and pitch motion is limited. Thus, the estimation of the vertical bias component is not sufficiently accurate to calibrate the magnetic field measurement.
In this paper, we propose that only the horizontal component of the magnetic field measurement should be used to prevent the adverse effect of the magnetic field measurement on the estimation of the roll and pitch, when the magnetic field is significantly distorted or the vertical bias component cannot be calibrated, owing to the limited roll and pitch motion.
The horizontal component of the magnetic field is the magnetic field projected onto the North-East plane (NE plane) in the NED coordinate frame. The proposed method extracts the horizontal component h m(t) of the magnetic field measurement from the raw magnetic field measurement m(t) by subtracting the vertical component of the magnetic field using Equation (46):

t) a(t) · a(t) a(t).
The last term of Equation (46) represents the projection of the magnetic field m(t) onto vector a(t), which is assumed to be perpendicular to the NE plane, as assumed when the roll and pitch are calculated from the acceleration measurement. The horizontal component h m(t) of the magnetic field measurement m(t) is normalized to h m u (t), as expressed by Equation (47): In this approach, the measurement consists of the normalized acceleration and normalized horizontal magnetic field, as expressed by Equation (48): Because the NE-plane magnetic field measurement is used in the measurement update procedure of the EKF, the NE-plane geomagnetic field h m g should be used instead of the full geomagnetic field m g , which is given by Equation (24). The NE-plane geomagnetic field h m g is simply the NE component of m g , as expressed by Equation (49): The use of h m(t) and h m g instead of m(t) and m g completes the approach where only the horizontal component of the magnetic field measurement is used to remove the adverse effect originating from the use of the full magnetic field measurement in the estimation of roll and pitch.

Experiments and Discussion
The proposed method was evaluated and compared with other methods through experiments inside two test pools with different remotely operated vehicles (ROVs) and sensors. The compared methods were NECF, SR, EKF, CF, and the method installed into the used AHRS by the manufacturer (LORD MicroStrain, Williston, VT, USA). The internal method is a complementary filter (CF) and provides the attitude in Euler angles [26]. In the two experiments, the attitude and location were represented in the NED coordinate frame. The origin of the NED coordinate frame for the location description is the initial location, where the UUV starts a predefined motion through a given trajectory for each test. The method proposed in this paper is called the FM method, which stands for field measurement method.
The first experiment uses an AHRS (LORD Microstrain 3DM-GX4 25) [26], a fibre optic gyroscope (FOG, Advanced Navigation, Spatial Fog, Sydney, NSW, Australia) [27], and a DVL (Teledyne RD Instruments, Navigator Doppler Velocity Log, Poway, CA, USA) [28]. The attitude provided by the FOG is used as a reference with which the attitude estimated by the abovementioned methods are compared. Figure 1 shows the installation of the AHRS and FOG in water-proof housing. The performance of the FOG is shown in Table 1. Table 2 describes the performance of the AHRS. Figure 2 shows the UUV used for the experiment along with the test pool. The test pool is located in the city of Gwangju, Korea, and has a length and depth of 5.5 m, 5.0 m, respectively. According to the world magnetic model for 2015-2020 at the location of the test pool, the geomagnetic field was set to m g = (0.307198G, −0.041355G, 0.390659G) T [25]. The UUV was controlled remotely to navigate through the predefined trajectories in the pool. The UUV had five degrees of freedom in motion with six thrusters: four horizontal thrusters for surge, sway, and yaw motion, and two vertical thrusters for roll and heave motion. Figure 3 shows the trajectory through which the UUV navigated within the test pool. The UUV moved in a circular trajectory from its initial location. Every time it reached back to its initial location at the end of one circular motion, it turned around in place to follow the reverse direction and trace back to the circular trajectory. The travel distance and time were 264 m, and 34 min and 20 s, with an average speed of 0.21 m/s when the UUV was moving, excluding the rotation at the end of each circulation.  Table 2. Specification of AHRS used in experiment 1.

Parameter Value
Gyroscope bias instability 10 • /h Accuracy of CF attitude estimation (static) ±0.5 • Accuracy of CF attitude estimation (dynamic) ±2.0 •   Table 3 shows the average (Mn), average absolute error (MnA), standard deviation (Std), and maximum peak to peak (PtP) of the error. In Table 3, the AHRS CF indicates the CF method implemented internally in the AHRS, while CF represents the CF implemented in this study. Figure 5 depicts the error statistics graphically to make the comparison easier. The mean, average absolute error, standard deviation, and peak to peak of the error are shown. Although the roll and pitch estimation errors for all methods are comparable, a difference in the yaw estimation error is noticeable. Therefore, the discussion regarding the results is focused on the statistical comparison of the yaw estimation error. Figure 5 shows a graphical representation of the mean, average absolute error, standard deviation, and peak to peak of the error. As shown in Table 3 and Figure 5, the proposed FM method achieves the best performance in peak to peak, and the second best performance with regard to the mean and average absolute error. The standard deviation performance of the FM ranks third. Although the performance of SR in the mean error and average absolute error are better than the performance of the proposed method, the performance of the standard deviation and peak to peak of the error for SR are inferior in comparison with the same performance of the proposed FM method. In comparison with the other methods, it can be generally considered that the FM method outperforms the other methods.     The second experiment used an AHRS (LORD Microstrain 3DM-GX3 25, LORD MicroStrain, Williston, VT, USA) [29] and a DVL (LinkQuest NavQuest 600 Micro, San Diego, CA, USA) [30] installed onto the UUV. The test pool was located in the city of Daejeon, Korea, and the geomagnetic field was set to m g = (0.300684G, −0.042526G, 0.401987G) T [25], in accordance with the 2015-2020 world magnetic model with regard to the specific location. The test pool and the UUV are shown in Figure 6. The UUV navigated through the two different trajectories shown in Figure 7-namely, a circular trajectory and a rectangular trajectory. For each trajectory, once the robot circulated and reached its initial location, it turned back and traced the trajectory to the initial location by following the reverse direction. Then, it turned back and traced the previous path, and moved again toward the initial location. The robot repeated this back and forth navigation, and finally returned to its initial location such that the destination of the navigation was the robot's initial location. The total travel distance of the circular motion was 282.7 m, and the total distance of the rectangular motion was 100.8 m [22]. Because the second experiment did not use a high-end attitude sensor, the attitude estimation performance was analyzed indirectly by comparing the accuracy of the dead-reckoned localization, which depends on the estimated attitude. Figure 8 shows the trajectories dead-reckoned using the estimated attitude and velocity measured by the DVL. Because the DVL measures the velocity in the sensor coordinate frame, the velocity is transformed to the velocity in the NED coordinate frame by using the attitude estimated by the considered methods. The NED velocity was dead-reckoned to produce the trajectory in the NED coordinate frame. Therefore, the accuracy of the calculated trajectory in the NED coordinate frame partially indicates the performance of attitude estimation. The UUV began moving from the origin (0, 0, 0) of the coordinate system, and ended its motion at the same location after circulating the trajectory. The errors in the estimated attitude and measured velocity cause the calculated UUV location to drift with time. The drift of the calculated final destination from the true destination (0, 0, 0) is regarded as a measure of accuracy for the estimated attitude. Table 4 describes the distance error of the estimated location at the final destination. Although the drift in the z-direction is not shown in Figure 8, the z-directional error adds to the distance error on the xy plane shown in Figure 8. Thus, it results in the distance error presented in Table 4. The distance error for the FM is shorter than the errors of the other methods. Although the distance error is not the exact measure of accuracy for the estimated attitude, the result proves that the calculated trajectory based on the proposed method is better than that of the other methods, and indirectly indicates that the proposed attitude estimation method is the most probable one amongst all methods compared in this study. Reduction of computational expenses are critical for the practical use of the algorithm since the computational capability is constrained by the limitations in space, weight, power consumption, and price available within the underwater vehicles. Table 5 compares the computation time for  the five methods. The table lists Table 5 indicates that the proposed method takes the longest computation time. However, all the algorithms run in comparable computation time, which is less than 0.20 ms. The proposed method can run over 5.0 KHz, and it is enough for real-time implementation. The NECF and EKFs have been widely used for underwater navigation. As the computation time for NECF and EKFs are acceptable for practical underwater navigations, it is expected that the computation time for the proposed method is also allowable for real-time application. Memory usage is another important factor for implementation of the algorithms for practical use. For the Matlab implementation, the proposed method uses 3864 bytes of variables, which is considered reasonable because the other methods also require comparable amount of variables with the proposed method.
In order to use the method for navigation in open water environment, the following aspects should be taken into account besides the computation expenses: (1) periodic position and attitude fixing is needed; and (2) geomagnetic field reference should be updated as the vehicle navigates from one place to another. Underwater navigation which depends on DVL and inertia measurement unit (IMU) inevitably suffer from accumulation of error in location and heading as shown in Figure 8. Therefore, periodic correction of the location and attitude is required. In some applications, the vehicle periodically surfaces and corrects the location and heading using the GPS data. In some other cases where the vehicle navigates in a limited navigation space, acoustic baseline system such as USBL, short baseline (SBL), or long baseline (LBL) can be used for the correction.
The proposed method requires geomagnetic field reference data, which depends on geographic location of the underwater vehicle. Table 6 shows the geomagnetic field values of two locations which are 9.1 km apart from each other. The geomagnetic field values are for the date of 9 January 2019, which are calculated using the world magnetic model (WMM) 2015, version 2, updated in September 2018 [31]. The distance between the two places are calculated using the Universal Transverse Mercator (UTM) coordinates with world geodetic system (WGS)-84 ellipsoid data. Table 6 indicates that the variation of geomagnetic field values is small enough to be ignored for application of the proposed method when the locations are within 9.1 km range. Table 6. Comparison of geomagnetic field of two locations 9.1 km apart, for the date of 9 January 2019. The geomagnetic field could be considered constant if the vehicle's workspace is within several kilometers of bound. However, if the vehicle navigates around significantly large range, the geomagnetic field should be updated. For the update, the geomagnetic model should be installed into the vehicle control processor system. Alternatively, the geomagnetic field value can be updated through wireless communication when the vehicle surfaces. Because the estimation error of location by the proposed method accumulates with time and travel distance, the geomagnetic field calculated based on the estimated location also deviates from the true value. Therefore, periodic fix of location is needed using either GPS reception at the surfacing period or using acoustic baseline systems.

Conclusions
This paper proposes a method of estimating the attitude of an underwater robot using measurement by an AHRS. The proposed method uses a magnetic field and acceleration measurement at the correction stage of the filtering. Unlike existing EKF based methods, the proposed method uses the field measurement as is, by directly comparing the measured field with the known field strength specific to the geographic locations where the UUV operates. The proposed method was compared with five other methods using measurement data sampled through two experiments in test tanks. The statistical analysis for the attitude error revealed that the proposed method performed better than, or at least comparably, with the other considered methods.
The proposed formulation uses Euler angles for the state transition and measurement model. It is expected that the use of an alternative representation for the attitude, such as a quaternion or rotation matrix, will result in a different state transition model and measurement model. In future work, an invariant observer which uses field measurement will be adopted [20] to further develop this study. It is thought that the incorporation of the field measurement approach with an alternative attitude representation utilizing the invariance property could improve the convergence performance and accuracy of attitude estimation.