# An Error Evaluation Method for the Ship Angular Flexure Measurement Based on the Principle of Relevance

Authors: Yutong Zhang; Xianglu Ma; Shiqiao Qin; Wei Wu; Wenfeng Tan
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Due to the lack of true ship angular flexure data, it is difficult to evaluate its measurement error of the angular velocity matching method in practice. In this paper, the cause of the measurement error of the ship flexure angle is analyzed in theory, and an evaluation method for the ship angular flexure measurement error based on the principle of relevance is proposed. The proposed method provides a prediction formula to describe the estimation error of the static flexure angle based on the off-diagonal elements of the error covariance matrixP in Kalman filtering. In addition, the optimized coefficient F is introduced to make the prediction error range better describe the real error variation. The optimized coefficient F ensures that the proposed formula has good prediction effects in all three directions. Simulations based on the actual measured ship flexure data are carried out, and the simulation results verify the capability of the prediction formula. The proposed method can be used in the evaluation of the ship flexure measurement error.


---

## Body

## 1. Introduction

Modern battleships are equipped with different kinds of sensors and weapon systems, such as radars, guided missiles, optoelectronic collimator, and other peripheral apparatus. The coordinate system misalignment caused by the ship flexure affects the performance of these systems. A unified spatial coordinate reference is urgently needed to transfer high-precision attitude information to these sensors and weapon systems.According to the researches of Day and Arruda [1] and Petovello et al. [2], the ship flexure measurement is the key procedure to establish a unified coordinate reference. Schneider [3] points out that large dynamic flexure may occur between the master inertial navigation system (MINS) and slave inertial navigation system (SINS) when the ship undergoes wave- or wind-induced loads or turning maneuvers due to the elastic structure of the ship. Moreover, there will be slow-varying flexure caused by the change of the ship’s mass distribution and temperature variation, which is named as the static flexure. The ship flexure measurement methods are widely studied by researchers all over the world, and the main solutions include the optical autocollimation method [4], photogrammetric method [5], GPS-based method [2], and the inertia vector matching method [6]. Among them, the inertia vector matching method can measure the static and dynamic flexures simultaneously without using external information, which makes it widely used in ship flexure measurement field.The inertia vector matching method uses two sets of inertial measurement units (IMUs) installed at different positions of the hull to provide observation information [3, 6]. The IMUs can measure the angular velocity vector and the acceleration vector of the ship at different positions. By constructing the inertia vector matching model, the ship flexure can be estimated with a Kalman filter [7]. In order to improve the accuracy of the inertia vector matching method, many researchers started to study the error mechanism of the ship flexure measurement. Zheng et al. point out that the ship flexure measurement error mainly consists of the method error and the instrumental error [8]. The instrumental error is caused by nonideal instruments such as the gyroscope bias. According to Zheng et al.’s study [9], this error is usually observable and can be compensated. However, the method error is irrelevant to the accuracy of instruments. It is caused by the incompleteness of the measurement equation. According to Wu et al.’s study [10], the angular velocity matching (AVM) equation is a kind of ill-conditioned equation without an analytical solution. The approach to determine the flexure angle is to take successive measurements and then estimate it using a Kalman filter. For the AVM method, the dynamic flexure component can be estimated accurately, while the static component is biased severely. The TA procedure has a large static estimation error even when the MINS and SINS are all equipped with high-quality gyro instruments and the dynamic flexure model parameters are determined. Therefore, an inherent measurement error exists. Many research studies [11–15] point out that this estimation error is caused by the coupling effect between the ship angular motion and the dynamic flexure. Browne et al. [11] suggested that the alignment error and estimation speed has a strong correlation with the ship angular motion and ship dynamic flexure. Wu et al. [12] explained the physical mechanism of the correlation. The dynamic flexure and angular motion are all the response of the ship structure to the wave loads, and they are likely to be correlated, which introduced a coupling error on the estimation of static flexure. However, the specific relationship between this coupling effect and the estimation error of static flexure is still unknown.Wu et al. [10] deduced a coupling error function based on the spatial geometric modelling and mathematical analysis. The function suggests that the coupling error value depends on the phase delay and amplitude ratio of the dynamic flexure and ship angular velocity. But the coupling error can only be estimated based on ship structural and hydrodynamic analysis, which is not given by the study. Therefore, the coupling error function mentioned in the study [10] is only a qualitative analysis and cannot be applied in the practical navigation process.In this paper, the estimation error of the static flexure angle, which is caused by the coupling between the ship angular motion and the dynamic flexure, is analyzed in theory. On the other hand, due to this coupling effect, it is difficult to measure the ship flexure accurately by the inertia vector matching method, and there is no clear report to evaluate its measurement accuracy so far. We propose an error evaluation method for the ship angular flexure measurement based on the principle of relevance. In order to better describe the error range, an optimized coefficientF is introduced to the proposed method. Simulations with real measured flexure data from a ship are implemented, and the results demonstrate the effectiveness of the proposed method.


## 2. Theory of AVM Method for Ship Angular Flexure Measurement

### 2.1. Definition of Coordinate Systems

The coordinate systems used in this paper are defined as follows.The nonrotating inertial coordinate system (i-frame) is fixed in inertial space and centered on the Earth. Its x-axis is in the equatorial plane and points to the vernal equinox. The z-axis is aligned with the Earth’s rotation axis and vertical to the equatorial plane. The y-axis completes a right-handed frame [16].The ship body coordinate system (b-frame) is rigidly attached to the ship and its origin Ob is at the ship’s mass center. The x-axis points to the starboard, y-axis is directed in the direction of ship bow, and z-axis points upward.The MINS coordinate system (m-frame) and the SINS coordinate system (s-frame) are attached to the MINS and SINS, respectively.


### 2.2. Angular Velocity Matching Equation

As shown in Figure1, the MINS and SINS, both of which contain laser gyroscope units (LGUs), are installed at different positions of the ship. If the ship is regarded as an absolutely rigid body, the angular velocities measured by the MINS and SINS with respect to i-frame satisfy the following equation:(1)Csmω→iss=ω→imm,where Csm denotes the transformation matrix from s-frame to m-frame, ω→imm denotes the angular velocity of m-frame with respect to i-frame, and ω→iss denotes the angular velocity of s-frame. Since the initial installation misalignment angle is usually small, it can be approximated as(2)Csm=I+φ→×,where I denotes the identity matrix with size 3 × 3 and φ→× is the skew-symmetric form of initial installation misalignment angle φ→, which can be expressed as(3)φ→×=0−φ→zφ→yφ→z0−φ→x−φ→yφ→x0.Figure 1
Schematic diagram of the measurement of the flexure angle.Substitute equation (2) into equation (1), and we have the following expression:(4)ω→iss−ω→imm=ω→iss×φ→.In fact, the ship is not an absolutely rigid body but has elastic flexure, which makes the misalignment Euler angle change with time. Therefore, equation (1) should be rewritten as(5)Csmω→iss=ω→imm+ω→msm,where ω→msm denotes the angular velocity from m-frame to s-frame. Since the flexure angle φ→ represents the Euler angle from m-frame to s-frame, φ→ and ω→msm satisfy the following equation:(6)φ→˙=ω→msm,where φ→˙ represents the angular flexure velocity and the dot operator •˙ represents the derivation with respect to time t. By substituting equations (5) and (6) into equation (4), we can get the AVM equation under elastic condition:(7)Δω→=ω→iss−ω→imm=ω→iss×φ→+φ→˙.Normally, the total misalignment angle is divided into two parts as described in the introduction part, which are named as the static componentφ→0 and the dynamic component φ→d, respectively:(8)φ→=φ→0+φ→d.Here the static componentφ→0 contains both the initial installation misalignment angle and slow-varying flexure angle. By substituting equation (8) into equation (7), we can get the following expression:(9)Δω→=ω→iss−ω→imm=ω→iss×φ→+φ→˙d.In fact, the gyros are not ideal. Considering the gyros bias, equation (9) can be expressed as(10)Δω→=ω→iss−ω→imm=ω→iss×φ→+φ→˙d+ε⟶,where ε⟶ is the gyro bias difference between MINS and SINS.


### 2.3. Optimal Estimation of Ship Flexure Angle

In Section “Angular velocity matching equation,” we establish the liner model ofΔω→ with respect to the flexure angle. The flexure angle can be estimated by a Kalman filter through establishing the state-space model and the observation model.The state vector is specified by(11)x=φ0xφ0yφ0zφdxφdyφdzφ˙dxφ˙dyφ˙dzεxεyεzT,where φ0xφ0yφ0zT, φdxφdyφdzT, and φ˙dxφ˙dyφ˙dzT are the subcomponents of φ→, φ→0, and φ→˙d, respectively. εxεyεzT is the vector of gyro bias difference between MINS and SINS. The dynamic component φ→d is generally modeled by the second-order Markov process [3, 6, 17, 18], which can be written as(12)φ→¨dj+2αjφ→dj+ω0j2φ→dj=2ω0jσjαjejt,j=x,y,z,where αj denotes the damping factor, σj denotes the standard deviation of φ→dj, ejt denotes Gaussian white noise with unit variance, and the circular frequency ω0j is given by(13)ω0j=αj2+βj2,where βj is the prevailing variation frequency of φ→dj.The state-space model can be expressed as(14)x˙=Fx+w,where the state matrix F is given by(15)F=03×303×303×303×303×303×3I3×303×303×3F32F3303×303×303×303×3I3×3,here, F32 and F33 satisfy(16)F32=−ω0x2000−ω0y2000−ω0z2,F33=−2αx000−2αy000−2αz.The 12 × 1 state noise vectorw can be expressed as(17)EwwT=diag0,…,0︸64ω0x2σx02αx4ω0y2σy2αy4ω0z2σz2αzσεx2σεy2σεz2,where E⋅ is the expectation operator.The measurement model can be expressed as(18)z=Hx+v,where z is the observation vector determined by equation (10), v is the measurement error vector, and H is the observation matrix. According to equation (10), H can be expressed as(19)H=ω→iss×ω→iss×03×3I3×3.Hence, the state-space model and the observation model are established. The state-space vector can be estimated by a Kalman filter and the optimal estimations of the flexure angle can be acquired.


## 2.1. Definition of Coordinate Systems

The coordinate systems used in this paper are defined as follows.The nonrotating inertial coordinate system (i-frame) is fixed in inertial space and centered on the Earth. Its x-axis is in the equatorial plane and points to the vernal equinox. The z-axis is aligned with the Earth’s rotation axis and vertical to the equatorial plane. The y-axis completes a right-handed frame [16].The ship body coordinate system (b-frame) is rigidly attached to the ship and its origin Ob is at the ship’s mass center. The x-axis points to the starboard, y-axis is directed in the direction of ship bow, and z-axis points upward.The MINS coordinate system (m-frame) and the SINS coordinate system (s-frame) are attached to the MINS and SINS, respectively.


## 2.2. Angular Velocity Matching Equation

As shown in Figure1, the MINS and SINS, both of which contain laser gyroscope units (LGUs), are installed at different positions of the ship. If the ship is regarded as an absolutely rigid body, the angular velocities measured by the MINS and SINS with respect to i-frame satisfy the following equation:(1)Csmω→iss=ω→imm,where Csm denotes the transformation matrix from s-frame to m-frame, ω→imm denotes the angular velocity of m-frame with respect to i-frame, and ω→iss denotes the angular velocity of s-frame. Since the initial installation misalignment angle is usually small, it can be approximated as(2)Csm=I+φ→×,where I denotes the identity matrix with size 3 × 3 and φ→× is the skew-symmetric form of initial installation misalignment angle φ→, which can be expressed as(3)φ→×=0−φ→zφ→yφ→z0−φ→x−φ→yφ→x0.Figure 1
Schematic diagram of the measurement of the flexure angle.Substitute equation (2) into equation (1), and we have the following expression:(4)ω→iss−ω→imm=ω→iss×φ→.In fact, the ship is not an absolutely rigid body but has elastic flexure, which makes the misalignment Euler angle change with time. Therefore, equation (1) should be rewritten as(5)Csmω→iss=ω→imm+ω→msm,where ω→msm denotes the angular velocity from m-frame to s-frame. Since the flexure angle φ→ represents the Euler angle from m-frame to s-frame, φ→ and ω→msm satisfy the following equation:(6)φ→˙=ω→msm,where φ→˙ represents the angular flexure velocity and the dot operator •˙ represents the derivation with respect to time t. By substituting equations (5) and (6) into equation (4), we can get the AVM equation under elastic condition:(7)Δω→=ω→iss−ω→imm=ω→iss×φ→+φ→˙.Normally, the total misalignment angle is divided into two parts as described in the introduction part, which are named as the static componentφ→0 and the dynamic component φ→d, respectively:(8)φ→=φ→0+φ→d.Here the static componentφ→0 contains both the initial installation misalignment angle and slow-varying flexure angle. By substituting equation (8) into equation (7), we can get the following expression:(9)Δω→=ω→iss−ω→imm=ω→iss×φ→+φ→˙d.In fact, the gyros are not ideal. Considering the gyros bias, equation (9) can be expressed as(10)Δω→=ω→iss−ω→imm=ω→iss×φ→+φ→˙d+ε⟶,where ε⟶ is the gyro bias difference between MINS and SINS.


## 2.3. Optimal Estimation of Ship Flexure Angle

In Section “Angular velocity matching equation,” we establish the liner model ofΔω→ with respect to the flexure angle. The flexure angle can be estimated by a Kalman filter through establishing the state-space model and the observation model.The state vector is specified by(11)x=φ0xφ0yφ0zφdxφdyφdzφ˙dxφ˙dyφ˙dzεxεyεzT,where φ0xφ0yφ0zT, φdxφdyφdzT, and φ˙dxφ˙dyφ˙dzT are the subcomponents of φ→, φ→0, and φ→˙d, respectively. εxεyεzT is the vector of gyro bias difference between MINS and SINS. The dynamic component φ→d is generally modeled by the second-order Markov process [3, 6, 17, 18], which can be written as(12)φ→¨dj+2αjφ→dj+ω0j2φ→dj=2ω0jσjαjejt,j=x,y,z,where αj denotes the damping factor, σj denotes the standard deviation of φ→dj, ejt denotes Gaussian white noise with unit variance, and the circular frequency ω0j is given by(13)ω0j=αj2+βj2,where βj is the prevailing variation frequency of φ→dj.The state-space model can be expressed as(14)x˙=Fx+w,where the state matrix F is given by(15)F=03×303×303×303×303×303×3I3×303×303×3F32F3303×303×303×303×3I3×3,here, F32 and F33 satisfy(16)F32=−ω0x2000−ω0y2000−ω0z2,F33=−2αx000−2αy000−2αz.The 12 × 1 state noise vectorw can be expressed as(17)EwwT=diag0,…,0︸64ω0x2σx02αx4ω0y2σy2αy4ω0z2σz2αzσεx2σεy2σεz2,where E⋅ is the expectation operator.The measurement model can be expressed as(18)z=Hx+v,where z is the observation vector determined by equation (10), v is the measurement error vector, and H is the observation matrix. According to equation (10), H can be expressed as(19)H=ω→iss×ω→iss×03×3I3×3.Hence, the state-space model and the observation model are established. The state-space vector can be estimated by a Kalman filter and the optimal estimations of the flexure angle can be acquired.


## 3. Coupling Effect between Ship Angular Motion and Dynamic Flexure

Due to the complexity of the ship’s mechanical model, it is difficult to analyze the coupling effect between the ship angular motion and dynamic flexure physically. This paper gives a statistical perspective to analyze this coupling effect, and the estimation error of the static flexure angle caused by this coupling effect is deduced in theory. And we further proposed an error evaluation method for the ship angular flexure measurement based on the principle of relevance.


### 3.1. Estimation Error of Static Flexure Angle Caused by Coupling Effect

According to the works of Kailath and Sorenson [19, 20, 21], the Kalman filter can be rightfully regarded as an efficient computational solution of the least square method. So it is reasonable to study the mechanism of bias error based on the ordinary least square (OLS) theory:(20)ΔΩ=Ωφ0+Θ˙,where(21)ΔΩ=Δω1Δω2⋮Δωk3k×1,Ω=ω1×ω2×⋮ωk×3k×3,Θ˙=φ˙d1φ˙d2⋮φ˙dk3k×1,here, k represents the ordinal number of sampling and Δω→× is the skew-symmetric form of vector Δω→. Let Θ˙ denote the real value of dynamic flexure velocity, Θ˙^ denote the estimation value of dynamic flexure angular velocity in Kalman filtering, and δΘ˙ denote the estimation error, which satisfy the following equation:(22)Θ˙=Θ˙^+δΘ˙.Substituting equation (22) into equation (20) results in the following expression:(23)ΔΩ−Θ˙^=Ωφ0+δΘ˙.If there is no estimation error of dynamic angular velocity (i.e.,δΘ˙=0), the static component of the flexure angle can be estimated by the OLS method with(24)φ^0=ΩTΩ−1ΩTΔΩ−Θ˙^.In fact,δΘ˙ cannot be ignored due to the estimation error of dynamic angular velocity. Substitute equation (23) into equation (24), and the least square estimation of φ^0 is as follows:(25)φ^0=ΩTΩ−1ΩTΩφ0+δΘ˙=ΩTΩ−1ΩTΩφ0+ΩTΩ−1ΩTδΘ˙=φ0+ΩTΩ−1ΩTδΘ˙.Therefore, the estimation error of the static flexure angle is given by(26)δφ0=φ0−φ^=−ΩTΩ−1ΩTδΘ˙.Using the sampling sequence of angular velocity data, equation (26) can be expressed as(27)δφ0=−ω1×Tω2×T⋯ωk×Tω1×Tω2×T⋯ωk×T−1ω1×Tω2×T⋯ω2×Tδφ˙d1δφ˙d2⋯δφ˙dk.Equation (27) can be expressed with the ordinal number of sampling k as(28)δφ0=−∑i=1kωy2i+ωz2ik−∑i=1kωxiωyik−∑i=1kωxiωzik−∑i=1kωyiωxik∑i=1kωx2i+ωz2ik−∑i=1kωyiωzik−∑i=1kωziωxik−∑i=1kωziωyik∑i=1kωx2i+ωy2ik−1∑i=1kωziδφ˙dyi−ωyiδφ˙dzik∑i=1kωxiδφ˙dzi−ωziδφ˙dxik∑i=1kωyiδφ˙dxi−ωxiδφ˙dyik.According to the Wiener–Khinchin law of large numbers,δφ0 will converge when k⟶∞, and the convergence value is given by(29)δφ0=A−1B=−Rωyωy0+Rωzωz0−Rωxωy0−Rωxωz0−Rωyωx0Rωxωx0+Rωzωz0−Rωyωz0−Rωzωx0−Rωzωy0Rωxωx0+Rωyωy0−1Rωzδφ˙dy0−Rωyδφ˙dz0Rωxδφ˙dz0−Rωzδφ˙dx0Rωyδφ˙dx0−Rωxδφ˙dy0.Equation (29) shows the estimation error of the static flexure angle caused by the coupling between angular velocity and dynamic flexure angular velocity. R represents the correlation function, Rωyωy0 represents the autocorrelation function of ωy without time delay, and Rωzδφ˙dy0 represents the cross-correlation function of ωz and δφ˙dy without time delay. The rest elements are of similar physical meanings. In equation (29), the diagonal elements of matrix A represent the autocorrelation functions of the ship’s angular velocities in the same direction while the nondiagonal elements represent cross-correlation functions in different directions. Since the value of autocorrelation function is much larger than the cross-correlation function for the same group of signals, equation (29) can be approximately simplified with diagonal elements reserved only and the subcomponents of the estimation error of the static flexure angle can be expressed as(30)δφ→0x=Rωyδφ˙dz0−Rωzδφ˙dy0Rωyωy0+Rωzωz0,δφ→0y=Rωzδφ˙dx0−Rωxδφ˙dz0Rωxωx0+Rωzωz0,δφ→0z=Rωxδφ˙dy0−Rωyδφ˙dx0Rωxωx0+Rωyωy0.In equation (30), the denominators are the autocorrelation function components, which can be calculated directly with the output data of IMUs. The numerators are the cross-correlation function between the ship angular velocity and the error of dynamic flexure angular velocity. Since the components of δφ→˙d in x, y, and z directions are unknown, equation (30) cannot be directly used to predict the estimation error of the flexure angle.


### 3.2. Error Evaluation Formula Based on the Principle of Relevance

According to the Kalman filter theory [22], the diagonal elements of covariance matrix P denote the variance of the state variable’s estimation error. And the estimation of state variable x^ subjects to Gauss distribution, of which x is the real value and P is the variance:(31)x^−x∼N0,P.According to the “three-sigma (3σ) rule” (Pauta criterion) in statistical theory, as P is the variance of x^−x, P can be regarded as 1σ range of the estimation error x^−x and 3P can be regarded as 3σ range. Because of the low observability of the static flexure angle (especially in z direction), it is unreasonable to describe the error range of δφ→0 with P^φ0,φ0 directly. On the other hand, the observability of dynamic flexure angular velocity is higher than that of the static flexure angle in z direction, so the element P^φ˙d,φ˙d can reflect the estimation error of δφ→˙d. However, the diagonal elements P^φ˙d,φ˙d are slow-varying elements and their values are always positive. If P^φ˙d,φ˙d is used to represent δφ→˙d in cross-correlation function, the prediction value in equation (30) will tend to zero. So the diagonal elements P^φ˙d,φ˙d cannot be used to represent δφ→˙d in the prediction formula. The off-diagonal elements in covariance matrix P denote the coupling error between two different state variables. Since the inherent measurement error of the static flexure angle is induced by the coupling between the ship angular velocity and dynamic flexure velocity, it is reasonable to evaluate the error of dynamic flexure angular velocity by off-diagonal element P^φo,φ˙d. Therefore, the estimation error of the static flexure angle in the Kalman filter can be derived by δφ→˙d with the following expression:(32)δφ→˙d,F=F⋅sgnPφ→o,φ→˙dPφ→o,φ→˙d,where F is the optimized coefficient. Normally, we set F = 3 to represent the 3σ range. P is the covariance matrix in the Kalman filter. sgn is the sign function, which is used to take the positive or negative sign of the off-diagonal element P^φo,φ˙d.Equation (32) shows the prediction error of the angular flexure velocity. Substitute equation (32) into equation (30), and the prediction error of the static flexure angle in three directions can be deduced as(33)δφ→0x,F=Rωyδφ˙dz,F0−Rωzδφ˙dy,F0Rωyωy0+Rωzωz0,δφ→0y,F=Rωzδφ˙dx,F0−Rωxδφ˙dz,F0Rωxωx0+Rωzωz0,δφ→0z,F=Rωxδφ˙dy,F0−Rωyδφ˙dx,F0Rωxωx0+Rωyωy0.Since the observability of dynamic flexure velocity is better than that of flexure angle, equation (33) has the potential to represent the error range with higher performance.


## 3.1. Estimation Error of Static Flexure Angle Caused by Coupling Effect

According to the works of Kailath and Sorenson [19, 20, 21], the Kalman filter can be rightfully regarded as an efficient computational solution of the least square method. So it is reasonable to study the mechanism of bias error based on the ordinary least square (OLS) theory:(20)ΔΩ=Ωφ0+Θ˙,where(21)ΔΩ=Δω1Δω2⋮Δωk3k×1,Ω=ω1×ω2×⋮ωk×3k×3,Θ˙=φ˙d1φ˙d2⋮φ˙dk3k×1,here, k represents the ordinal number of sampling and Δω→× is the skew-symmetric form of vector Δω→. Let Θ˙ denote the real value of dynamic flexure velocity, Θ˙^ denote the estimation value of dynamic flexure angular velocity in Kalman filtering, and δΘ˙ denote the estimation error, which satisfy the following equation:(22)Θ˙=Θ˙^+δΘ˙.Substituting equation (22) into equation (20) results in the following expression:(23)ΔΩ−Θ˙^=Ωφ0+δΘ˙.If there is no estimation error of dynamic angular velocity (i.e.,δΘ˙=0), the static component of the flexure angle can be estimated by the OLS method with(24)φ^0=ΩTΩ−1ΩTΔΩ−Θ˙^.In fact,δΘ˙ cannot be ignored due to the estimation error of dynamic angular velocity. Substitute equation (23) into equation (24), and the least square estimation of φ^0 is as follows:(25)φ^0=ΩTΩ−1ΩTΩφ0+δΘ˙=ΩTΩ−1ΩTΩφ0+ΩTΩ−1ΩTδΘ˙=φ0+ΩTΩ−1ΩTδΘ˙.Therefore, the estimation error of the static flexure angle is given by(26)δφ0=φ0−φ^=−ΩTΩ−1ΩTδΘ˙.Using the sampling sequence of angular velocity data, equation (26) can be expressed as(27)δφ0=−ω1×Tω2×T⋯ωk×Tω1×Tω2×T⋯ωk×T−1ω1×Tω2×T⋯ω2×Tδφ˙d1δφ˙d2⋯δφ˙dk.Equation (27) can be expressed with the ordinal number of sampling k as(28)δφ0=−∑i=1kωy2i+ωz2ik−∑i=1kωxiωyik−∑i=1kωxiωzik−∑i=1kωyiωxik∑i=1kωx2i+ωz2ik−∑i=1kωyiωzik−∑i=1kωziωxik−∑i=1kωziωyik∑i=1kωx2i+ωy2ik−1∑i=1kωziδφ˙dyi−ωyiδφ˙dzik∑i=1kωxiδφ˙dzi−ωziδφ˙dxik∑i=1kωyiδφ˙dxi−ωxiδφ˙dyik.According to the Wiener–Khinchin law of large numbers,δφ0 will converge when k⟶∞, and the convergence value is given by(29)δφ0=A−1B=−Rωyωy0+Rωzωz0−Rωxωy0−Rωxωz0−Rωyωx0Rωxωx0+Rωzωz0−Rωyωz0−Rωzωx0−Rωzωy0Rωxωx0+Rωyωy0−1Rωzδφ˙dy0−Rωyδφ˙dz0Rωxδφ˙dz0−Rωzδφ˙dx0Rωyδφ˙dx0−Rωxδφ˙dy0.Equation (29) shows the estimation error of the static flexure angle caused by the coupling between angular velocity and dynamic flexure angular velocity. R represents the correlation function, Rωyωy0 represents the autocorrelation function of ωy without time delay, and Rωzδφ˙dy0 represents the cross-correlation function of ωz and δφ˙dy without time delay. The rest elements are of similar physical meanings. In equation (29), the diagonal elements of matrix A represent the autocorrelation functions of the ship’s angular velocities in the same direction while the nondiagonal elements represent cross-correlation functions in different directions. Since the value of autocorrelation function is much larger than the cross-correlation function for the same group of signals, equation (29) can be approximately simplified with diagonal elements reserved only and the subcomponents of the estimation error of the static flexure angle can be expressed as(30)δφ→0x=Rωyδφ˙dz0−Rωzδφ˙dy0Rωyωy0+Rωzωz0,δφ→0y=Rωzδφ˙dx0−Rωxδφ˙dz0Rωxωx0+Rωzωz0,δφ→0z=Rωxδφ˙dy0−Rωyδφ˙dx0Rωxωx0+Rωyωy0.In equation (30), the denominators are the autocorrelation function components, which can be calculated directly with the output data of IMUs. The numerators are the cross-correlation function between the ship angular velocity and the error of dynamic flexure angular velocity. Since the components of δφ→˙d in x, y, and z directions are unknown, equation (30) cannot be directly used to predict the estimation error of the flexure angle.


## 3.2. Error Evaluation Formula Based on the Principle of Relevance

According to the Kalman filter theory [22], the diagonal elements of covariance matrix P denote the variance of the state variable’s estimation error. And the estimation of state variable x^ subjects to Gauss distribution, of which x is the real value and P is the variance:(31)x^−x∼N0,P.According to the “three-sigma (3σ) rule” (Pauta criterion) in statistical theory, as P is the variance of x^−x, P can be regarded as 1σ range of the estimation error x^−x and 3P can be regarded as 3σ range. Because of the low observability of the static flexure angle (especially in z direction), it is unreasonable to describe the error range of δφ→0 with P^φ0,φ0 directly. On the other hand, the observability of dynamic flexure angular velocity is higher than that of the static flexure angle in z direction, so the element P^φ˙d,φ˙d can reflect the estimation error of δφ→˙d. However, the diagonal elements P^φ˙d,φ˙d are slow-varying elements and their values are always positive. If P^φ˙d,φ˙d is used to represent δφ→˙d in cross-correlation function, the prediction value in equation (30) will tend to zero. So the diagonal elements P^φ˙d,φ˙d cannot be used to represent δφ→˙d in the prediction formula. The off-diagonal elements in covariance matrix P denote the coupling error between two different state variables. Since the inherent measurement error of the static flexure angle is induced by the coupling between the ship angular velocity and dynamic flexure velocity, it is reasonable to evaluate the error of dynamic flexure angular velocity by off-diagonal element P^φo,φ˙d. Therefore, the estimation error of the static flexure angle in the Kalman filter can be derived by δφ→˙d with the following expression:(32)δφ→˙d,F=F⋅sgnPφ→o,φ→˙dPφ→o,φ→˙d,where F is the optimized coefficient. Normally, we set F = 3 to represent the 3σ range. P is the covariance matrix in the Kalman filter. sgn is the sign function, which is used to take the positive or negative sign of the off-diagonal element P^φo,φ˙d.Equation (32) shows the prediction error of the angular flexure velocity. Substitute equation (32) into equation (30), and the prediction error of the static flexure angle in three directions can be deduced as(33)δφ→0x,F=Rωyδφ˙dz,F0−Rωzδφ˙dy,F0Rωyωy0+Rωzωz0,δφ→0y,F=Rωzδφ˙dx,F0−Rωxδφ˙dz,F0Rωxωx0+Rωzωz0,δφ→0z,F=Rωxδφ˙dy,F0−Rωyδφ˙dx,F0Rωxωx0+Rωyωy0.Since the observability of dynamic flexure velocity is better than that of flexure angle, equation (33) has the potential to represent the error range with higher performance.


## 4. Simulations and Results

### 4.1. Verification of the Proposed Method

The 19-hour measured ship flexure data of “Yuanwang No. 3” [23] with the large steel pipe-based method, which can be regarded as the real flexure data, are utilized to validate the proposed method. The angular velocities measured by the MINS and SINS are from actual navigation systems installed on the ship. In the simulation, the static flexure angle is set to [360, 360, 360] The 19-hour dynamic ship flexure data are shown in Figure 2, and the dynamic flexure angle can be treated as Markov processes with parameters listed in Table 1. The performance of gyros is shown in Table 2.Figure 2
Actual ship flexure data used for simulation.Table 1
Model parameters of dynamic flexure angle.
σ (arcs)
ω
0
α (1/s)
X, Y, Z-axes
0.0033
0.4189
0.1Table 2
The performance of gyros.
Bias°/h
Random walk°/h
X, Y, Z-axes
0.003
0.001The outputs of MINS and SINS are then used to compute the ship flexure. The estimation results are shown in Figure3. The blue curves represent the real values of the flexure angle in three directions, and the red curves represent the estimation results.Figure 3
Estimation results of ship flexure in different directions: (a)x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)The estimation error of the flexure angle is shown with green curves in Figure4, the 3σ error range of δφ→0 calculated with P^φ0,φ0 directly is shown with blue curves in Figure 4(a), and the 3σ error range with the proposed method δφ→0,F is shown with red curves in Figure 4(b). Obviously, the 3σ error range calculated with P^φ0,φ0 cannot conclude the real estimation error. However, the estimation error fluctuates within the 3σ error range calculated with the proposed method, which coincides with the analysis in Section “Error evaluation formula based on the principle of relevance.”Figure 4
Prediction 3σ error range of δφ→0 described by (a) P^φ0,φ0 and (b) δφ→0,F.
(a)
(b)As analyzed in Section “Error evaluation formula based on the principle of relevance,” since the observability of dynamic flexure angular velocity is higher than that of the static flexure angle, the 3σ range of the proposed equation (33) has better performance to describe the estimation error. When F = 3, the components of dynamic flexure angular velocity in x and y directions are shown in Figure 5. The blue curves represent the real values of δφ→˙dx and δφ→˙dy, respectively, and the red curves represent the prediction range. Obviously, the curves of the prediction error range in both x and y directions are consistent with the fluctuation phase of the real estimation error, and the amplitudes are slightly greater than that of the real estimation error, which coincides with the analysis in Section “Error evaluation formula based on the principle of relevance.”Figure 5
Prediction results of dynamic flexure angular velocity in (a)x direction and (b) y direction.
(a)
(b)


### 4.2. Analysis of the Optimized CoefficientF

In order to further study the influence ofF on the prediction performance and get the most appropriate F, we analyze the prediction performance with different values of F. Figure 6 shows the relationship between the prediction error range of the static flexure angle in x, y, and z directions and the value of F, respectively. The symbols on the red line are the sampling data points, which represent the mean values of 19-hour results. It can be seen that the prediction error range increases linearly as F getting greater in each direction. Therefore, by adjusting the value of F, we can control the prediction error range to approach to the real estimation error.Figure 6
The relationship between the prediction error range andF: (a) x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)Figure7 shows the relationship between the prediction error range and the real estimation error in x, y, and z directions when F changes from 0 to 2. The real estimation error is only related to time t and independent of F. However, the prediction error range is related to both time t and F. The intersection of these two surfaces defines a curve, on which the prediction error range is equal to the real estimation error.Figure 7
The surfaces of the prediction error range and the real estimation error: (a)x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)Normally, the prediction error range should be greater than (or equal to) the real estimation error, so as to ensure that the real error value falls within the estimation range. At the same time, the prediction error range should not be too large for the purpose of avoiding overestimation. Therefore, the meaningful prediction error surface is the part that is little higher than the intersection curve. According to our simulations, the most appropriate value of optimized coefficientF is [1.4, 0.8, 2]T, and the prediction error range is shown in Figure 8. It can be seen that the curve of real estimation error fluctuates within the prediction error range in every direction. The magnitude of real estimation error reaches the maximum at time epoch nearing the 12 h, which intersects the curve of the prediction error range. Therefore, the measurement error falls within the prediction range, which strongly validates the effectiveness of the proposed method.Figure 8
Prediction results of the static flexure angle withF = [1.4, 0.8, 2]T: (a) x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)


## 4.1. Verification of the Proposed Method

The 19-hour measured ship flexure data of “Yuanwang No. 3” [23] with the large steel pipe-based method, which can be regarded as the real flexure data, are utilized to validate the proposed method. The angular velocities measured by the MINS and SINS are from actual navigation systems installed on the ship. In the simulation, the static flexure angle is set to [360, 360, 360] The 19-hour dynamic ship flexure data are shown in Figure 2, and the dynamic flexure angle can be treated as Markov processes with parameters listed in Table 1. The performance of gyros is shown in Table 2.Figure 2
Actual ship flexure data used for simulation.Table 1
Model parameters of dynamic flexure angle.
σ (arcs)
ω
0
α (1/s)
X, Y, Z-axes
0.0033
0.4189
0.1Table 2
The performance of gyros.
Bias°/h
Random walk°/h
X, Y, Z-axes
0.003
0.001The outputs of MINS and SINS are then used to compute the ship flexure. The estimation results are shown in Figure3. The blue curves represent the real values of the flexure angle in three directions, and the red curves represent the estimation results.Figure 3
Estimation results of ship flexure in different directions: (a)x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)The estimation error of the flexure angle is shown with green curves in Figure4, the 3σ error range of δφ→0 calculated with P^φ0,φ0 directly is shown with blue curves in Figure 4(a), and the 3σ error range with the proposed method δφ→0,F is shown with red curves in Figure 4(b). Obviously, the 3σ error range calculated with P^φ0,φ0 cannot conclude the real estimation error. However, the estimation error fluctuates within the 3σ error range calculated with the proposed method, which coincides with the analysis in Section “Error evaluation formula based on the principle of relevance.”Figure 4
Prediction 3σ error range of δφ→0 described by (a) P^φ0,φ0 and (b) δφ→0,F.
(a)
(b)As analyzed in Section “Error evaluation formula based on the principle of relevance,” since the observability of dynamic flexure angular velocity is higher than that of the static flexure angle, the 3σ range of the proposed equation (33) has better performance to describe the estimation error. When F = 3, the components of dynamic flexure angular velocity in x and y directions are shown in Figure 5. The blue curves represent the real values of δφ→˙dx and δφ→˙dy, respectively, and the red curves represent the prediction range. Obviously, the curves of the prediction error range in both x and y directions are consistent with the fluctuation phase of the real estimation error, and the amplitudes are slightly greater than that of the real estimation error, which coincides with the analysis in Section “Error evaluation formula based on the principle of relevance.”Figure 5
Prediction results of dynamic flexure angular velocity in (a)x direction and (b) y direction.
(a)
(b)


## 4.2. Analysis of the Optimized CoefficientF

In order to further study the influence ofF on the prediction performance and get the most appropriate F, we analyze the prediction performance with different values of F. Figure 6 shows the relationship between the prediction error range of the static flexure angle in x, y, and z directions and the value of F, respectively. The symbols on the red line are the sampling data points, which represent the mean values of 19-hour results. It can be seen that the prediction error range increases linearly as F getting greater in each direction. Therefore, by adjusting the value of F, we can control the prediction error range to approach to the real estimation error.Figure 6
The relationship between the prediction error range andF: (a) x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)Figure7 shows the relationship between the prediction error range and the real estimation error in x, y, and z directions when F changes from 0 to 2. The real estimation error is only related to time t and independent of F. However, the prediction error range is related to both time t and F. The intersection of these two surfaces defines a curve, on which the prediction error range is equal to the real estimation error.Figure 7
The surfaces of the prediction error range and the real estimation error: (a)x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)Normally, the prediction error range should be greater than (or equal to) the real estimation error, so as to ensure that the real error value falls within the estimation range. At the same time, the prediction error range should not be too large for the purpose of avoiding overestimation. Therefore, the meaningful prediction error surface is the part that is little higher than the intersection curve. According to our simulations, the most appropriate value of optimized coefficientF is [1.4, 0.8, 2]T, and the prediction error range is shown in Figure 8. It can be seen that the curve of real estimation error fluctuates within the prediction error range in every direction. The magnitude of real estimation error reaches the maximum at time epoch nearing the 12 h, which intersects the curve of the prediction error range. Therefore, the measurement error falls within the prediction range, which strongly validates the effectiveness of the proposed method.Figure 8
Prediction results of the static flexure angle withF = [1.4, 0.8, 2]T: (a) x direction; (b) y direction; (c) z direction.
(a)
(b)
(c)


## 5. Conclusions

In this paper, the cause of the measurement error of the ship flexure angle is analyzed in theory, and it is proved that the cross-coupling between the ship angular motion and the dynamic flexure leads to the estimation error of the static flexure angle. We further propose an evaluation method for the measurement error of the ship flexure angle based on the covariance matrix of Kalman filtering. According to the covariance matrix of Kalman filtering of the AVM method, the prediction error range formula of the static flexure angle error is derived. The optimized coefficientF is introduced to make the prediction error range better describe the real estimation error. Simulations with actual ship flexure data show the effectiveness of the proposed prediction formula of flexure angle estimation error range. And the prediction formula has the best performance with F = 1.6 in x direction, F = 1 in y direction, and F = 2 in z direction. Therefore, the optimized coefficient F can be set as a vector to ensure that the prediction error range formula has good prediction effects in all three directions. The proposed method can be used in the evaluation of the ship flexure measurement error.



---

*Source: 9795650-2019-11-11.xml*