Particle Filter with Novel Nonlinear Error Model for Miniature Gyroscope-Based Measurement While Drilling Navigation

The derivation of a conventional error model for the miniature gyroscope-based measurement while drilling (MGWD) system is based on the assumption that the errors of attitude are small enough so that the direction cosine matrix (DCM) can be approximated or simplified by the errors of small-angle attitude. However, the simplification of the DCM would introduce errors to the navigation solutions of the MGWD system if the initial alignment cannot provide precise attitude, especially for the low-cost microelectromechanical system (MEMS) sensors operated in harsh multilateral horizontal downhole drilling environments. This paper proposes a novel nonlinear error model (NNEM) by the introduction of the error of DCM, and the NNEM can reduce the propagated errors under large-angle attitude error conditions. The zero velocity and zero position are the reference points and the innovations in the states estimation of particle filter (PF) and Kalman filter (KF). The experimental results illustrate that the performance of PF is better than KF and the PF with NNEM can effectively restrain the errors of system states, especially for the azimuth, velocity, and height in the quasi-stationary condition.


Introduction
Horizontal drilling (HD), especially reentry multilateral horizontal well drilling (RMHWD), receives considerable attention because it provides great economic value in the oil and gas industry.The RMHWD technology was originally proposed in an industry consortium of operators and service companies in 1997 [1] and RMHWD can significantly increase the production of oil and gas by extending the contact areas between oil and gas reservoirs and the well drilling equipment and by revisiting the existing wellbores [2].RMHWD is also readily available under challenging drilling sites such as offshore zones, mountain areas, and downtown centers [3].The branches of the drill pipe in RMHWD radiate from the main wellbore to the targeted reservoir sections [4].The high pressure waterjet drilling technology promotes the development of RMHWD since the conventional large rotary drill bit can be replaced by a small-diameter waterjet drill bit.The conventional measurement-while-drilling (MWD) tool combines three-axis gyroscopes and three-axis magnetometers to provide the attitude (tool face, inclination, and azimuth) of the drill bit.However, the performance of the MWD tools would deteriorate under the strong magnetic interference caused by the following factors: (1) flow of currents in the atmosphere and moving of solar wind [5]; (2) existing drilling fluid and debris [6,7]; (3) ferromagnetic deposits in the proximity of bottom hole assembly (BHA) [8].Non-magnetic drill collars can effectively eliminate the magnetic disturbance; however, the high cost, high weight, and weak fracture drill collars fail to be widely used.
Alternatively, gyroscope-based MWD (GWD) tools can continuously provide the attitude, velocity, and position of the drill bit regardless of the strong magnetic interferences.GWD tools with pure strapdown inertial navigation algorithm were proposed in [9,10] and the tools can save considerable rig time and remove the blind during the kickoff [11].The size of the fiber optic gyroscope or the laser optic gyroscope limits their application to the surveying of downhole drilling.The emergence of a low-cost, small-size, and high-performance MEMS sensor brings the gyroscopes to the small-diameter RMHWD.The MGWD device with two-axis gyroscopes and three-axis accelerometers in the application of small-diameter space (less than 24 millimeters) was proposed in [12], in which the angular rate (Y-axis) of MGWD tools is calculated by the accelerometers in the horizontal plane when the drilling equipment is under quasi-stationary conditions.The MGWD system can provide high-precision navigation solutions in short-term [13] and can perform real-time data transmission between the MGWD tools and the ground operation center.However, the integration calculation within resolving the mechanization of navigation solutions accumulates large errors because of biases of gyroscopes and accelerometers.The errors of inertial sensors are composed of two parts: deterministic errors and stochastic errors [14].The deterministic errors can partially be reduced by the calibration approach.Therefore, the residual of deterministic errors and the stochastic errors are main error sources of inertial sensors.To sum up, the problems should be solved for MGWD system are the elimination of the residual of deterministic errors and stochastic errors.
The stochastic errors of inertial sensors can be modeled by the first-order Gaussian-Markov (GM) model or the high-order autoregressive (AR) model [15].The Levenberg-Marquardt iterative least squares can fit the nonlinear parametric model of stochastic errors as well [16].The traditional error models of the MGWD system, such as the Phi-angle or Psi angle model, are based on the assumption that the attitude errors are small enough that the DCM can be simplified by the small-angle attitude errors.The linear modified error model for large heading uncertainty was proposed in [17] and the model of moving alignment with large errors was put forward in [18].In the application of the Kalman filter (KF), both error models neglect the high-order error terms that may introduce errors when the attitude errors are large.In addition, the system noises and the observation noises are not Gaussian for MGWD system in sophisticated downhole drilling environments (high temperature and high pressure).Therefore, the KF is suboptimal compared with the particle filter (PF), which is the recursive approximation of the posterior probability density function (pdf) with particles and weights [19].The conventional in-drilling alignment (IDA) approaches include the zero velocity update (ZUPT) alignment, velocity matching alignment, rotary modulation alignment, etc.The ZUPT alignment can limit the position errors up to 40 m during a 90-min experiment by using a Litton LTN90-100 inertial measurement unit (IMU) [20].The velocity matching alignment merges velocity from global positioning system (GPS) or other velocity sensors to the filtering for compensating the errors [21].An IDA approach was introduced to restrict the azimuth error of the drill bit in [22,23]; the azimuth error was reduced 25 times smaller than the conventional magnetometer based MWD tools.The rotary-in-drilling alignment (RIDA) was proposed in terms of the theory of rotary modulation alignment in an inertial navigation system (INS) to minimize the dynamic position errors and the azimuth error of the MWD system [24,25].However, GPS signals cannot reach to a deep downhole and the size of IDA or RIDA is beyond the dimension requirements of the RMHWD.
This paper proposes a NNEM for the MGWD system under large-angle attitude error conditions.The error matrix of DCM is introduced to represent the attitude errors.Because of the nonlinearity of the system state space model (SSM), the PF is selected as the state estimation approach.The body frame is the right, forward, and up direction of the MGWD system and the navigation frame is the north, east, and down direction in this paper.
In Section 2, the NNEM of the MGWD system is proposed.Sections 3 and 4 introduce the recursive Bayesian estimation theory and PF algorithm, respectively.Results and experiments are explained in Section 5. Eventually, conclusions are presented in Section 6.

NNEM of the MGWD System
The error model of the strapdown inertial navigation system (SINS), based on the discrete time differential equation and the error model, plays a fundamental role in the estimation of system states.Conventional linear error models (LEM) of SINS, such as the Phi-angle model, the Psi-angle model [26,27], and the quaternion error model [28] are widely used in aerospace and marine navigation.The Phi-angle error model analyzes the SINS in the true navigation frame while the Psi-angle model is performed in the computer frame.In this paper, the Psi-angle model is selected for the states estimation of KF while the error model in [29] is used as the NEM for the states estimation of the PF.
The derivation of a linear model by the perturbation method alone can be implemented in light of the assumption that the attitude errors are small enough.The initial alignment of the MGWD system cannot provide an accurate attitude because of the low performance of the MEMS sensors and the sophisticated downhole drilling environments.Therefore, this paper proposes an accurate NNEM of the MGWD system under large-angle attitude error conditions.
The key variable for the derivation of the SINS error model is the rotation vector, which is the eigenvector of the DCM associated with eigenvalue 1.The eigenvector of the DCM can be expressed as follows [30,31]: where R n b is the DCM from the body frame to the navigation frame, Φ is the rotation vector, and I is a 3 ˆ3 identity matrix.The rotation vector is proposed in [32], and the relationship between rotation vector and DCM can be expressed as: where pΦˆq is the skew symmetric matrix of Φ, i.e.,: where Φ x , Φ y , and Φ z are the components of Φ projected in X-axis, Y-axis, and Z-axis.The magnitude of Φ equals: So far, the problem of seeking an attitude solution equates to seeking a solution for the rotation vector.The derivation of the kinematic equation of rotation vector was reviewed in [33], and the conventional differential equation of rotation vector was given by [34]: where ω b nb is the angular rate of the body frame with respect to the navigation frame resolved in the body frame.
This section reviews the nonlinear error model of MGWD according to [27], in which the attitude error model is derived according to the linear quaternion error model.The attitude errors are defined in terms of the rotation vector: Φ " pδφ, δθ, δψq where δφ, δθ, and δψ are the roll error, pitch error, and azimuth error.The NNEM of the attitude can be expressed as follows: .
where c n and s n are defined parameters that can be expressed as: where ω b ib is the angular rate obtained from the gyroscope, ω n in is the angular rate of the navigation frame with respect to the inertial frame resolved in the navigation frame, and ω n in can be expressed as: where ω n ie is the Earth rotation rate in the navigation frame and ω n ie can be written as follows: where ω ie is the Earth rotation rate 7.292115 ˆ10 ´5rad{s « 15.041067 ˝{h; ϕ is the latitude; ω n en is the angular rate of the navigation frame with respect to the Earth frame resolved in the navigation frame; and ω n en can be written as: where V E and V N represent the east velocity and north velocity, respectively, in the navigation frame; h is the ellipsoidal height; R N and R M are the radii of curvature in the meridian and the prime vertical, respectively; and R N , R M can be expressed as follows [35]: where R , e are the semi-major axis and linear eccentricity of the reference ellipsoid, respectively, in which R " 6378137 m and e " 0.0818191908426.The output errors of accelerometers and gyroscopes (δf b ib and δω b ib ) can be expressed as [36]: where ∆SF x,g and ∆SF z,g are the scale factor errors of X-axis gyroscope and Z-axis gyroscope.∆SF y,g is defined as the scale factor error of a virtual gyroscope along the Y-axis.∆δ xy,g , ∆δ xz,g , ∆δ yz,g , ∆δ yx,g , ∆δ zy,g , and ∆δ zx,g are non-orthogonal gyroscope errors along the X-axis, Y-axis, and Z-axis.
∆SF x,a , ∆SF y,a , and ∆SF z,a are the scale factor errors of the X-axis accelerometer, Y-axis accelerometer, and Z-axis accelerometer, respectively.∆δ xy,a , ∆δ xz,a , ∆δ yz,a , ∆δ yx,a , ∆δ zy,a , and ∆δ zx,g are non-orthogonal acceleration errors along the X-axis, Y-axis, and Z-axis.x,g , y,g , z,g and x,a , y,a , z,a are the angular drift errors of gyroscopes and drift errors of accelerometers along the X-axis, Y-axis, and Z-axis, respectively.
Through differentiating the Equations ( 11) and ( 12), δω n in can be written as [37]: where δV E and δV N are the errors of east velocity and north velocity, respectively; δϕ and δh are the errors of latitude and height, respectively; and δω n ie and δω n en are the errors of ω n ie and ω n en , respectively.The conventional velocity update equation can be written as follows [21]: where V n is the velocity in the navigation frame; f b ib is the output of accelerometers; and g n is the gravity expressed in the navigation frame (the model of g n is proposed in [38]).
The velocity error model can be obtained by perturbing Equation (20), and the DCM can be approximated by the following expression: where the designator " in the above variables represents variables with errors.Ξ is a skew symmetric matrix of attitude error and Ξ is given by: However, the above equation may introduce errors if large attitude errors exist, so the error matrix of DCM p∆R n b q is introduced for the derivation of velocity error equation.The DCM in Equation ( 21) can be replaced by: The velocity update equation with error interferences can be written as: .
V n ¯can be obtained by the difference between Equations ( 20) and (24), ignoring the second-order error product R b n ∆R n b δf b ib and `2δω n ie `δω n en ˘δV n , yields: The only unknown parameter in Equation ( 25) is ∆R n b , which can be obtained in terms of Φ, and the relationship between ∆R n b and Φ can be expressed as [27]: ´s2 n ´pδθq 2 `pδψq 2 ¯sn p1 `cn q δψ `s2 n δφδθ ´sn p1 `cn q δθ `s2 n δφδψ ´sn p1 `cn q δψ `s2 n δφδθ ´s2 n ´pδφq 2 `pδψq 2 ¯sn p1 `cn q δφ `s2 n δθδψ s n p1 `cn q δθ `s2 n δφδψ ´sn p1 `cn q δφ `s2 The velocity error equation is obtained by inserting Equation (26) into Equation (25).There is no relationship between the model position error ´δ .P ¯and the attitude, so the position error model is still linear as [39]: where δV D is the error of the down velocity.
The nonlinear error model of the MGWD system can be achieved through Equations ( 7), ( 25) and (27).The observation equation is the linear ZUPT equation, which is expressed as follows: where Z k is the difference between the estimate state pV N , V E , V D , ϕ, λ, hq and the zero velocity and zero position (V N,0 , V E,0 , V D,0 , ϕ 0 , λ 0 , h 0 ): H " where λ is the longitude and V N,0 , V E,0 , V D,0 , ϕ 0 , λ 0 , and h 0 represent the zero velocity and zero position.H 12 , H 13 , H 21 , and H 23 are a 3 ˆ3 zero matrix; H 11 , H 22 are a 3 ˆ3 identity matrix.

Recursive Bayesian Estimation
In this section, we briefly review the theory of the recursive Bayesian estimation.The system is generally modeled as an SSM, which includes a system equation and an observation equation in the discrete time domain; the SSM can be written as follows: x k " f k´1 px k´1 , k´1 q (31) where f k´1 : R n ˆRm Ñ R n is the system nonlinear function; k´1 P R m is a zero mean, white noise with known pdf p p k´1 q; x k and x k´1 are the system state vectors at present time k and previous moment k-1; z k P R n is the measurement from external reference; h k : R n ˆRr Ñ R p is the measurement function; and v k P R r is a zero mean, white noise with known pdf p pv k q. Figure 1 illustrates the recursive Bayesian estimation model in graphical form [40,41].

Recursive Bayesian Estimation
In this section, we briefly review the theory of the recursive Bayesian estimation.The system is generally modeled as an SSM, which includes a system equation and an observation equation in the discrete time domain; the SSM can be written as follows: where f k−1 : ℝ n × ℝ m → ℝ n is the system nonlinear function; ϖ k−1 ∈ ℝ m is a zero mean, white noise with known pdf p(ϖ k−1 ); x k and x k−1 are the system state vectors at present time k and previous moment k-1; z k ∈ ℝ n is the measurement from external reference; h k : ℝ n × ℝ r → ℝ p is the measurement function; and v k ∈ ℝ r is a zero mean, white noise with known pdf p(v k ). Figure 1 illustrates the recursive Bayesian estimation model in graphical form [40,41].A more general system model (Markov model) and observation model corresponding to the Equations ( 31) and ( 32) can be given by the following expressions: The definition of SSM in Equations ( 33) and ( 34) are based on the prior function and likelihood function.Both models are equivalent, but the general model is more convenient for the theoretical derivation.
The assumptions of recursive Bayesian estimation are that (1) the state vector x k is a hidden Markov process with initial distribution p(x 0 ) and x k is independent of the other state vectors Recursive Bayesian estimation model.The first horizontal row represents the states of the system, which propagate in terms of Equation (31), and the column represents the measurement process, which is formulated by Equation (32).x 0 , x 1 , x k´2 , x k´1 , x k are the system states from time 0 to time k and z 1 , z k´2 , z k´1 , z k are the measurements from time 1 to time k.
A more general system model (Markov model) and observation model corresponding to the Equations ( 31) and ( 32) can be given by the following expressions: The definition of SSM in Equations (33) and (34) are based on the prior function and likelihood function.Both models are equivalent, but the general model is more convenient for the theoretical derivation.
The assumptions of recursive Bayesian estimation are that (1) the state vector x k is a hidden Markov process with initial distribution p px 0 q and x k is independent of the other state vectors except x k´1 , i.e., the conditional pdf p px k |x 0:k´1 q " p px k |x k´1 q; (2) the measurement process z k only depends upon x k , i.e., the conditional pdf p pz k |x k , z 1:k´1 q " p pz k |x k q.
To sum up, the evaluation of the Bayesian estimation is given by the following steps.
(1) Prediction: evaluating the prior pdf p `xk |z 1:k´1 ˘by the knowledge of the observation before time k in the light of the Chapman-Kolmogorov equation [42]: where the pdf of initial state ppx 0 |z 0 q can be assigned by p px 0 |z 0 q " p px 0 q.p `xk |x k´1 ˘is the prior distribution, which is given by [43]: where δ p¨q is the Dirac delta function.
(2) Updating: the posterior pdf p px k |z 1:k q can be updated with the new measurement z k at time k in terms of the Bayes rule as follows [44]: Sensors 2016, 16, 371 8 of 24 where p pz k |x k q is the likelihood function, which is written by the following equation [43]: (3) State estimation: once the posterior is obtained, the estimation of state px k q and the covariance matrix of estimation error (P xk ) are given by [45]: xk where r x k " x k ´x k is the estimation error and E p¨q is the expectation of the random variables.

PF Algorithm
Because the function p `zk |z 1:k´1 ˘in Equation (37) cannot typically be computed analytically, sequential Monte Carlo (SMC) approximation, the key idea of the PF, was proposed for the hidden states estimation of the dynamic systems [46,47].The PF uses a finite set of weighted particles or samples to reconstruct the posterior pdf p px 0:k |z 1:k q of the state vectors regardless of the nonlinear or non-Gaussian systems [48,49].The implementation of the PF can be summarized in three main steps: sequential importance sampling (SIS), resampling, and roughening.The SIS plays a fundamental rule in the iteration of the PF since searching proper importance function decides the success of PF.Resampling aims at alleviating the particles' or samples' degeneracy and the roughing approach can be utilized to avoid the particles' or samples' impoverishment.

SIS Algorithm
The SIS algorithm steps from the literature [50,51].The basic idea of SIS is the selection of particles or samples with weights to represent the complicated high dimensional distributions [52].The posterior pdf ppx 0:k |z 1:k q is approximated by the random particles or samples and weights ) at time k as: where N is the number of particles and x 0:k is the system state.
The first task to implement the SIS is the selection of the importance function π px 0:k |z 1:k q , which is proportional of the posterior pdf ppx 0:k |z 1:k q and is easy to draw samples compared with ppx 0:k |z 1:k q.The weights can be updated in terms of the following equation: Generally, the importance function is selected such that [53]: Note that the Equation ( 41) can be factored as following equation according to the Bayesian rule [54]: Therefore, the weights can be updated by substituting Equations ( 43) and ( 44) into Equation ( 42), yielding: where x i k is the variable obtained by the evaluation of x i k´1 in the system model and x i k´1 is the samples drawn from the importance function π ´xi k |x i k´1 , z k ¯, which is related to the previous state x i k´1 and current observation z k .In the application of the MGWD system, the estimation of the current state is more important than the previous states, so Equation ( 41) can be rewritten as: As indicated in Equations ( 45) and ( 46), the posterior pdf depends on the likelihood function, the priori pdf, and the importance function.Since the likelihood function and the priori pdf can be obtained in terms of the SSM, the choice of importance function is essential for PF.The optimal function, the likelihood function, and the priori function are the components of the importance function [54].
The optimal function p ´xi k |x i k´1 , z k ¯was proposed as the importance function in [41]: Substituting Equation (47) into Equation ( 45) yields: The drawbacks of using optimal importance function are: (1) it is difficult to draw samples from the optimal function p ´xi k |x i k´1 , z k ¯; (2) it needs the integration calculation to update the weights.The optimal function can be factorized as the following format based on the Bayesian rule: As indicated in Equation ( 49), p ´xi k |x i k´1 , z k ¯9p `zk |x i k ˘if the likelihood function p `zk |x i k ˘is more noisy and it is integral in x i k [55].In this case, the importance function can be presented by the likelihood function, namely: Then, the weight in Equation ( 45) can be updated by the following equation: However, the likelihood function as the importance function can only be used in the high signal-to-noise ratio (SNR) case.In this paper, the priori function is selected as the importance function because of the tradeoff between computational complexity and the feasibility of implementation.The importance function is written by: Sensors 2016, 16, 371 10 of 24 Then, the weights can be updated by: The above equation shows that the weights only depend on the previous weights and the likelihood function, which are available from the knowledge of the observation equation and the pdf of the observation noise.In this application, we assume that the measurement noise of the MGWD system is additive, Gaussian white noise, so the likelihood function can be evaluated in a simple way [56]: where R is the covariance matrix of measurement noise v k .The normalization of the particle weights are performed to guarantee the sum of the particle weights equal to 1, and the normalization of the weights can be expressed as follows: Under the premise that the initial pdf ppx 0 |z 0 q is known, drawing particles x i 0 " ppx 0 |z 0 q, where i " 1, 2, N, as: where x i 0 is the initial particles and i 0 is the noise vector which is drawn from the pdf of system noise p p k´1 q.The initial weights are given by: The SIS algorithm can only be implemented in special cases because the variance of the importance weights only increase over time; eventually, they give rise to particle degeneracy.Particle degeneracy means that the updating particles are useless for the approximation of a posterior pdf.In order to solve the particles degeneracy problem, a resampling method is described in next section.

Resampling Algorithm
The PF would suffer from the particle degeneracy without a resampling algorithm.Figure 2 shows the principle of the particle resampling.
Sensors 2016, 16, 371 9 of 23 because of the tradeoff between computational complexity and the feasibility of implementation.The importance function is written by: Then, the weights can be updated by: The above equation shows that the weights only depend on the previous weights and the likelihood function, which are available from the knowledge of the observation equation and the pdf of the observation noise.In this application, we assume that the measurement noise of the MGWD system is additive, Gaussian white noise, so the likelihood function can be evaluated in a simple way [56]: where R is the covariance matrix of measurement noise v k .The normalization of the particle weights are performed to guarantee the sum of the particle weights equal to 1, and the normalization of the weights can be expressed as follows: Under the premise that the initial pdf p(x 0 |z 0 ) is known, drawing particles x 0 i ~p(x 0 |z 0 ), where i = 1, 2, N, as: where x 0 i is the initial particles and ϖ 0 i is the noise vector which is drawn from the pdf of system noise p(ϖ k−1 ).The initial weights are given by: The SIS algorithm can only be implemented in special cases because the variance of the importance weights only increase over time; eventually, they give rise to particle degeneracy.Particle degeneracy means that the updating particles are useless for the approximation of a posterior pdf.In order to solve the particles degeneracy problem, a resampling method is described in next section.

Resampling Algorithm
The PF would suffer from the particle degeneracy without a resampling algorithm.Figure 2 shows the principle of the particle resampling.The key of the resampling algorithm is that particles with large weights are replicated and particles with small weights are discarded [57].The effective sample size N eff is an indicator of the degree of the particles' degeneracy and a metric of employing the particles' resampling.N eff is provided by the following equation [58]: where var ´wpiq k ¯is the variation of the weight.The particles should be resampled if N eff is less than resample threshold pN th q, which equals 4N{5 in this paper.Equation (58) shows that N eff approaches 1 if the particles are severely degenerated; on the other hand, N eff approaches N if the particles have equal weight p1{Nq.The number of the resampled particles is not necessarily the same as the propagated particles N; however, the number of particles can be the same for computational convenience.Therefore, the total resampling particles are unchanged, still N, so the resampled particles have the same weight p1{Nq.
There are several ways to implement the resampling algorithm such as systematic resampling [59], multinomial resampling [60], residual resampling [61], etc.The resampling strategy in PF is sampling with replacement, that is, the particles drawn from the box will be replaced after the decimation.In this paper, the multinomial resampling approach is selected to alleviate particles' degeneracy.The multinomial resampling can be considered as a series of Bernoulli trials.N random samples u pnq , n " 1, 2, ¨¨¨N are generated from the uniform distribution over p0, 1s and the probabilities of resampling the particles is the same as u pnq when the following inequality is satisfied: Then, the nth particle x pnq k with the nth weight is the selection of the new samples or particles x piq k of the resampling.

Roughing Strategy
The problem of resampling is that only a few particles are drawn, which is also referred to as samples impoverishment.Therefore, the particles lack diversities because the selection may include many repeated particles.In order to address the samples impoverishment problem, the random noise is added to the particles after resampling.The posterior particles x piqk are written as: where ρ piq k is a zero mean random noise.At this point, the particles after the roughing step can be propagated for the update of the weights and the posterior pdf in the next time cycle.The first two moments of estimation derived from the particles can be approximated by the following equations:

Results and Discussion
The objectives of the designed experiments are to demonstrate that: (1) the performance of the PF with NEM or NNEM is better than KF with LEM under small-angle attitude error condition (Experiment 1); (2) PF with NNEM can provide more accurate states estimation for the MGWD system with large-angle attitude error (Experiment 2).The configuration of the MGWD device and the details of the experimental data collection are briefly described in the first part of Section 5.

MGWD Device for Validation of Experimental Results
The MGWD device consists of the azimuth MEMS sensor, the pitch MEMS sensor, and the microcontroller.Figure 3 shows the configuration of the MGWD device.Azimuth MEMS sensor provides the angular rate of the drill bit along the Z-axis while the angular rate along the X-axis is obtained from the pitch MEMS sensor.For implementation of the pure INS algorithm, the angular rate along the Y-axis is evaluated by the specific force obtained from the pitch MEMS sensor [12].

Results and Discussion
The objectives of the designed experiments are to demonstrate that: (1) the performance of the PF with NEM or NNEM is better than KF with LEM under small-angle attitude error condition (Experiment 1); (2) PF with NNEM can provide more accurate states estimation for the MGWD system with large-angle attitude error (Experiment 2).The configuration of the MGWD device and the details of the experimental data collection are briefly described in the first part of Section 5.

MGWD Device for Validation of Experimental Results
The MGWD device consists of the azimuth MEMS sensor, the pitch MEMS sensor, and the microcontroller.Figure 3 shows the configuration of the MGWD device.Azimuth MEMS sensor provides the angular rate of the drill bit along the Z-axis while the angular rate along the X-axis is obtained from the pitch MEMS sensor.For implementation of the pure INS algorithm, the angular rate along the Y-axis is evaluated by the specific force obtained from the pitch MEMS sensor [12].The MEMS sensor (SCC1300-D04, Murata Manufacturing Co., Ltd, Kyoto, Japan) is a combined gyroscope and accelerometer component and the sensor is characterized by low cost, small diameter, and high performance.The microcontroller (CC2538, Texas Instruments, Dallas, TX, USA) is a highspeed System-on-Chip (SoC) combined with ARM cortex-M3 processor.Tables 1 and 2 provide the technical specification of the MEMS sensor and the microcontroller, respectively.The MEMS sensor (SCC1300-D04, Murata Manufacturing Co., Ltd, Kyoto, Japan) is a combined gyroscope and accelerometer component and the sensor is characterized by low cost, small diameter, and high performance.The microcontroller (CC2538, Texas Instruments, Dallas, TX, USA) is a high-speed System-on-Chip (SoC) combined with ARM cortex-M3 processor.Tables 1 and 2 provide the technical specification of the MEMS sensor and the microcontroller, respectively.The types of data communication between the MEMS sensor and microcontroller is the serial peripheral interfaces (SPI) interface with 115,200 bit/s data transmission rate.The data collection software evaluates the navigation solutions with angular rate and specific force received from MGWD device through the universal serial bus (USB) with 20 Hz sample rate.
Figure 4 shows the process of the data collection and the contents of the experiments.The computer continuously reads the angular rate and the specific force of the drill bit from the static MGWD device.Therefore, the zero velocity and zero position are readily available as the reference of experiments, designed in terms of the size of attitude error.The types of data communication between the MEMS sensor and microcontroller is the serial peripheral interfaces (SPI) interface with 115,200 bit/s data transmission rate.The data collection software evaluates the navigation solutions with angular rate and specific force received from MGWD device through the universal serial bus (USB) with 20 Hz sample rate.
Figure 4 shows the process of the data collection and the contents of the experiments.The computer continuously reads the angular rate and the specific force of the drill bit from the static MGWD device.Therefore, the zero velocity and zero position are readily available as the reference of experiments, designed in terms of the size of attitude error.In the experiment, we collect original data from the static MGWD device for the period of 10 min to demonstrate the performance of the PF and the accuracy of the proposed NNEM. Figure 5 shows the original data from two-axis gyroscopes and three-axis accelerometers.The outputs of gyroscopes include the noise and the angular rate of the Earth projected onto the X-axis and Z-axis of the MGWD device.The X-axis and the Y-axis outputs of accelerometers are only the noise since the device is under static condition while the Z-axis output of accelerometer consists of the noise and gravity.In the experiment, we collect original data from the static MGWD device for the period of 10 min to demonstrate the performance of the PF and the accuracy of the proposed NNEM. Figure 5 shows the original data from two-axis gyroscopes and three-axis accelerometers.The outputs of gyroscopes include the noise and the angular rate of the Earth projected onto the X-axis and Z-axis of the MGWD device.The X-axis and the Y-axis outputs of accelerometers are only the noise since the device is under static condition while the Z-axis output of accelerometer consists of the noise and gravity.The types of data communication between the MEMS sensor and microcontroller is the serial peripheral interfaces (SPI) interface with 115,200 bit/s data transmission rate.The data collection software evaluates the navigation solutions with angular rate and specific force received from MGWD device through the universal serial bus (USB) with 20 Hz sample rate.
Figure 4 shows the process of the data collection and the contents of the experiments.The computer continuously reads the angular rate and the specific force of the drill bit from the static MGWD device.Therefore, the zero velocity and zero position are readily available as the reference of experiments, designed in terms of the size of attitude error.In the experiment, we collect original data from the static MGWD device for the period of 10 min to demonstrate the performance of the PF and the accuracy of the proposed NNEM. Figure 5 shows the original data from two-axis gyroscopes and three-axis accelerometers.The outputs of gyroscopes include the noise and the angular rate of the Earth projected onto the X-axis and Z-axis of the MGWD device.The X-axis and the Y-axis outputs of accelerometers are only the noise since the device is under static condition while the Z-axis output of accelerometer consists of the noise and gravity.

Experiment 1: Comparison of KF with LEM and PF with NEM under Small-Angle Attitude Error Condition
The aim of Experiment 1 is to confirm that the performance of the PF with NEM is better than the performance of KF with LEM under small-angle attitude error conditions.The assumption of KF is that all the states are Gaussian and the SSM is linear.However, the complicated drilling environments introduce such strong interference to the MGWD device that the linear SSM and Gaussian states are impossible.
The first 50 s are utilized for initial alignment; the results of the initial alignment for roll, pitch, and azimuth are ´0.022108˝, 0.180922 ˝, and ´10.351435 ˝, which are also set as the initial attitude in this experiment.Therefore, the initial attitude error is so small that the DCM can be approximated by the attitude.The rest of parameters related to the KF and PF are given as follows:

‚
Initial covariance matrix of system noise: R = diag([0.The aim of Experiment 1 is to confirm that the performance of the PF with NEM is better than the performance of KF with LEM under small-angle attitude error conditions.The assumption of KF is that all the states are Gaussian and the SSM is linear.However, the complicated drilling environments introduce such strong interference to the MGWD device that the linear SSM and Gaussian states are impossible. The first 50 s are utilized for initial alignment; the results of the initial alignment for roll, pitch, and azimuth are −0.022108°,0.180922°, and −10.351435°, which are also set as the initial attitude in this experiment.Therefore, the initial attitude error is so small that the DCM can be approximated by the attitude.The rest of parameters related to the KF and PF are given as follows:  Initial covariance matrix of system noise: R = diag([0.Figure 8 presents the evaluation of state density in terms of the updated pdf; the results illustrates that the pdf of states are partially Gaussian.Therefore, the KF is suboptimal since the assumption is that all states are Gaussian in KF while the PF would not be limited by the non-Gaussian assumption.Figure 8 presents the evaluation of state density in terms of the updated pdf; the results illustrates that the pdf of states are partially Gaussian.Therefore, the KF is suboptimal since the assumption is that all states are Gaussian in KF while the PF would not be limited by the non-Gaussian assumption.Figure 9 gives a performance comparison of the KF with LEM and PF with NEM or NNEM under small-angle attitude error.Figure 9 a1, a2 show that both the latitude and the longitude estimated by KF and PF converge to 45.7776° and 126.6879°, which are the same as the reference values (initial latitude and initial longitude).However, the height estimated by KF is 132.0782m while the height estimated by PF is 123.9999m, which is almost the same as the reference height (124 m).Therefore, the NEM or NNEM can restrain the errors of height effectively under small-angle attitude error conditions.(a1) (a2) (b2) (b2) Figure 9 gives a performance comparison of the KF with LEM and PF with NEM or NNEM under small-angle attitude error.Figure 9 a 1 , a 2 show that both the latitude and the longitude estimated by KF and PF converge to 45.7776 ˝and 126.6879 ˝, which are the same as the reference values (initial latitude and initial longitude).However, the height estimated by KF is 132.0782m while the height estimated by PF is 123.9999m, which is almost the same as the reference height (124 m).Therefore, the NEM or NNEM can restrain the errors of height effectively under small-angle attitude error conditions.Figure 9 gives a performance comparison of the KF with LEM and PF with NEM or NNEM under small-angle attitude error.Figure 9 a1, a2 show that both the latitude and the longitude estimated by KF and PF converge to 45.7776° and 126.6879°, which are the same as the reference values (initial latitude and initial longitude).However, the height estimated by KF is 132.0782m while the height estimated by PF is 123.9999m, which is almost the same as the reference height (124 m).Therefore, the NEM or NNEM can restrain the errors of height effectively under small-angle attitude error conditions.In Figure 9 b1, b2 the north velocity, east velocity, and down velocity estimated by PF are close to zero (0.0002 m/s, 0.0001 m/s, and 0.0001 m/s) while the errors of north velocity, east velocity, and down velocity estimated by KF are up to 0.3702 m/s, 0.2839 m/s, and 2.2953 m/s respectively.Therefore, the PF with NEM or NNEM gives better performance than KF with LEM does for the velocity estimation under small-angle attitude error.
In Figure 9 c1, c2 the roll, pitch, and azimuth of MGWD system estimated by PF with NEM or NNEM largely converge to −0.022125°, 0.180922°, and −10.353632°, respectively, which are the same as the reference points provided by the initial alignment, while the results of pitch and azimuth estimated by KF deviate from the reference values by about 8° and 10°, respectively.Therefore, the PF with NEM or NNEM can restrain the pitch and azimuth error effectively.

Experiment 2: Validation of the Performance of NNEM with Large-Angle Attitude Error
The objective of Experiment 2 is the validation of the performance of NNEM with large-angle attitude error.The experimental data used in this section are the same as those in Experiment 1 and the initial parameter settings are also the same as in experiment 1 apart from the initial attitude.The assumption of this experiment is that the initial alignment fails to provide an accurate initial attitude.We set the initial roll, pitch, and azimuth as 0°.The experiment is practical since the initial attitude is generally difficult to obtain in the application of downhole drilling.
Figure 10 a1, a2, b1, b2, c1, c2 show the estimation of roll, pitch, and azimuth of the MGWD system and the corresponding estimation errors.The roll, pitch, and azimuth converge to −0.021976°, 0.184392°, and −9.893302°, respectively, which are close to the reference values for initial alignment.The selection of initial attitude has little influence on the convergence of the attitude estimation by PF with NNEM since the roll converges to the reference value at around 320 s, the pitch converges to the reference value at around 180 s, and the azimuth converges to the reference value at around 360 s, as illustrated in Figure 10 a1, b1, and c1. Figure 10 a2, b2, c2 show that the error of roll and error of pitch estimated by PF with both NNEM and NEM are less than 0.0001° and 0.001°.The error of azimuth estimated by PF with NEM is 2.459532° while the estimation error of azimuth with NNEM is 2.5 × 10 −5 °.In addition, the accurate solutions of roll and pitch presented in Figure 10 demonstrate that the error of roll and error of pitch couple with the accuracy of east velocity and north velocity.Therefore, the PF with NNEM can greatly restrain the azimuth error under large-angle latitude error conditions.In Figure 9 b 1 , b 2 the north velocity, east velocity, and down velocity estimated by PF are close to zero (0.0002 m/s, 0.0001 m/s, and 0.0001 m/s) while the errors of north velocity, east velocity, and down velocity estimated by KF are up to 0.3702 m/s, 0.2839 m/s, and 2.2953 m/s respectively.Therefore, the PF with NEM or NNEM gives better performance than KF with LEM does for the velocity estimation under small-angle attitude error.
In Figure 9 c 1 , c 2 the roll, pitch, and azimuth of MGWD system estimated by PF with NEM or NNEM largely converge to ´0.022125 ˝, 0.180922 ˝, and ´10.353632 ˝, respectively, which are the same as the reference points provided by the initial alignment, while the results of pitch and azimuth estimated by KF deviate from the reference values by about 8 ˝and 10 ˝, respectively.Therefore, the PF with NEM or NNEM can restrain the pitch and azimuth error effectively.

Experiment 2: Validation of the Performance of NNEM with Large-Angle Attitude Error
The objective of Experiment 2 is the validation of the performance of NNEM with large-angle attitude error.The experimental data used in this section are the same as those in Experiment 1 and the initial parameter settings are also the same as in experiment 1 apart from the initial attitude.The assumption of this experiment is that the initial alignment fails to provide an accurate initial attitude.We set the initial roll, pitch, and azimuth as 0 ˝.The experiment is practical since the initial attitude is generally difficult to obtain in the application of downhole drilling.The error of azimuth estimated by PF with NEM is 2.459532 ˝while the estimation error of azimuth with NNEM is 2.5 ˆ10 ´5˝.In addition, the accurate solutions of roll and pitch presented in Figure 10 demonstrate that the error of roll and error of pitch couple with the accuracy of east velocity and north velocity.Therefore, the PF with NNEM can greatly restrain the azimuth error under large-angle latitude error conditions.Figure 11 a1, a2, b1, b2, c1, c2 illustrate the results of the north velocity, east velocity, and down velocity and the corresponding estimation errors of the MGWD system.The velocity estimated by PF, regardless of NEM or NNEM, apparently converges to zero velocity while the velocity error estimated by PF with NNEM is approximately 10 times smaller than KF with LEM.Therefore, the velocity errors estimated by PF with NNEM can be reduced effectively.

Conclusions
This paper proposes a NNEM for MGWD system under large-angle attitude error conditions.The error of DCM as the intermediate variable is introduced for the derivation of NNEM, which can effectively eliminate the approximation of errors existing in the derivation of conventional NEM.The evaluation of the state density shows that the system states are partially Gaussian, so the KF is

Conclusions
This paper proposes a NNEM for MGWD system under large-angle attitude error conditions.The error of DCM as the intermediate variable is introduced for the derivation of NNEM, which can effectively eliminate the approximation of errors existing in the derivation of conventional NEM.The evaluation of the state density shows that the system states are partially Gaussian, so the KF is suboptimal in this application.PF can improve the estimation performance of nonlinear problems by approximating the posterior pdf with the particles and weights.The PF with NEM gives a better performance than KF with LEM for the estimation of height, velocity, pitch, and azimuth in Experiment 1.The NNEM of the MGWD system provides an accurate estimation of position, velocity, and attitude compared with PF with NEM and KF with LEM in Experiment 2.

Figure 1 .
Figure 1.Recursive Bayesian estimation model.The first horizontal row represents the states of the system, which propagate in terms of Equation (31), and the column represents the measurement process, which is formulated by Equation (32).x 0 , x 1 , x k−2 , x k−1 , x k are the system states from time 0 to time k and z 1 , z k−2 , z k−1 , z k are the measurements from time 1 to time k.

1 Figure 1 .
Figure 1.Recursive Bayesian estimation model.The first horizontal row represents the states of the system, which propagate in terms of Equation (31), and the column represents the measurement process, which is formulated by Equation (32).x 0 , x 1 , x k´2 , x k´1 , x k are the system states from time 0 to time k and z 1 , z k´2 , z k´1 , z k are the measurements from time 1 to time k.

Figure 3 .
Figure 3. Configuration of the MGWD device.Each MEMS sensor is composed of a single-axis gyroscope and three-axis accelerometers.The body frame is the right (X-axis), forward (Y-axis) and up (Z-axis) direction.

Temperature − 40 Figure 3 .
Figure 3. Configuration of the MGWD device.Each MEMS sensor is composed of a single-axis gyroscope and three-axis accelerometers.The body frame is the right (X-axis), forward (Y-axis) and up (Z-axis) direction.

Figure 4 .
Figure 4. Process of data collection and contents of experiments.The MEMS sensors sense the angular rate and specific force of the drill bit and then transfer the data to the computer through a USB interface for evaluation of navigation solutions.The proposed algorithm would be verified by two groups of designed experiments with the reference of zero velocity and zero position.The design rule of experiments is based on the initial attitude error.

Figure 5 .
Figure 5. Original data obtained from experiments using the MGWD device.(a) Original data from two-axis gyroscopes; (b) original data from three-axis accelerometers.The length of experimental data is 600 s in quasi-stationary conditions.

Figure 4 .
Figure 4. Process of data collection and contents of experiments.The MEMS sensors sense the angular rate and specific force of the drill bit and then transfer the data to the computer through a USB interface for evaluation of navigation solutions.The proposed algorithm would be verified by two groups of designed experiments with the reference of zero velocity and zero position.The design rule of experiments is based on the initial attitude error.

Figure 4 .
Figure 4. Process of data collection and contents of experiments.The MEMS sensors sense the angular rate and specific force of the drill bit and then transfer the data to the computer through a USB interface for evaluation of navigation solutions.The proposed algorithm would be verified by two groups of designed experiments with the reference of zero velocity and zero position.The design rule of experiments is based on the initial attitude error.

Figure 5 .
Figure 5. Original data obtained from experiments using the MGWD device.(a) Original data from two-axis gyroscopes; (b) original data from three-axis accelerometers.The length of experimental data is 600 s in quasi-stationary conditions.

Figure 5 .
Figure 5. Original data obtained from experiments using the MGWD device.(a) Original data from two-axis gyroscopes; (b) original data from three-axis accelerometers.The length of experimental data is 600 s in quasi-stationary conditions.

Figure 6
Figure 6 illustrates the initial particles and initial weights.

Figure 6 .
Figure 6.Initial particles p(x 0 |z 0 ) and initial weights w 0 i for the PF.The initial particles and weights of attitude errors (δϕ, δθ, δψ), velocity errors (δV N , δV E , δV D ), and position errors (δφ, δλ, δh) are identical because the initial parameters (initial states and initial estimation errors) to evaluate the initial particles and initial weights are the same.

Figure 7
Figure 7 a1, a2, b1, b2, c1, c2 show the particles of the state errors together with the corresponding weights drawn from the priori pdf {p(x k i |x k−1 i )} at time k.The posterior pdf {p(x k |z k )} can be evaluated by the Equation (46) after the particles and the corresponding weights are obtained.The number of the particles and weights for each state are 200.The particles and weights are updated at every time circle; meanwhile, the posterior pdf is updated.

Figure 6 .
Figure 6.Initial particles p px 0 |z 0 q and initial weights w i 0 for the PF.The initial particles and weights of attitude errors pδφ, δθ, δψq, velocity errors pδV N , δV E , δV D q, and position errors pδϕ, δλ, δhq are identical because the initial parameters (initial states and initial estimation errors) to evaluate the initial particles and initial weights are the same.

Figure 7 a 1 ,Figure 7 .
Figure 7 a 1 , a 2 , b 1 , b 2 , c 1 , c 2 show the particles of the state errors together with the corresponding weights drawn from the priori pdf !p ´xi k |x i k´1 ¯) at time k.The posterior pdf tp px k |z k qu can be evaluated by the Equation (46) after the particles and the corresponding weights are obtained.The number of the particles and weights for each state are 200.The particles and weights are updated at every time circle; meanwhile, the posterior pdf is updated.

Figure 7 .
Figure 7. Particles or samples, along with the corresponding weights for evaluation of the posterior function.(a 1 ) and (a 2 ) are particles of velocity errors and the corresponding weights; (b 1 ) and (b 2 ) are the particles of attitude errors and the corresponding weights; (c 1 ) and (c 2 ) are the particles of positon errors and the corresponding weights.The priori pdf is the importance density function.The designator δL, the same as δϕ, is the latitude error.

Figure 8 .
Figure 8. Evaluation of state density.

Figure 8 .
Figure 8. Evaluation of state density.

Figure 8 .
Figure 8. Evaluation of state density.

Figure 9 .
Figure 9.Comparison of the performance of KF with LEM and PF with NEM or NNEM under small-angle attitude errors conditions.(a 1 ) and (a 2 ) are the estimated attitude and the corresponding estimation error; (b 1 ) and (b 2 ) are the estimated velocity and the corresponding estimation error; (c 1 ) and (c 2 ) are the estimated position and the corresponding estimation error.Plotting the figure uses the double Y-axis plotting approach because the results are in different orders of magnitude.

Figure 10 a 1 ,
a 2 , b 1 , b 2 , c 1 , c 2 show the estimation of roll, pitch, and azimuth of the MGWD system and the corresponding estimation errors.The roll, pitch, and azimuth converge to ´0.021976 ˝, 0.184392 ˝, and ´9.893302 ˝, respectively, which are close to the reference values for initial alignment.The selection of initial attitude has little influence on the convergence of the attitude estimation by PF with NNEM since the roll converges to the reference value at around 320 s, the pitch converges to the reference value at around 180 s, and the azimuth converges to the reference value at around 360 s, as illustrated in Figure 10 a 1 , b 1 , and c 1 .Figure 10 a 2 , b 2 , c 2 show that the error of roll and error of pitch estimated by PF with both NNEM and NEM are less than 0.0001 ˝and 0.001 ˝.

Figure 10 .
Figure 10.Results of attitude and estimation error of attitude by KF with LEM, PF with NEM, and PF with NNEM under large-angle attitude error conditions.(a1) and (a2) are the estimation of roll and estimation error of roll; (b1) and (b2) are the estimation of pitch and estimation error of pitch; (c1) and (c2) are the estimation of azimuth and estimation error of azimuth.Plotting the figure uses the triple Y-axis plotting approach because the results are in different orders of magnitude.

Figure 10 .
Figure 10.Results of attitude and estimation error of attitude by KF with LEM, PF with NEM, and PF with NNEM under large-angle attitude error conditions.(a 1 ) and (a 2 ) are the estimation of roll and estimation error of roll; (b 1 ) and (b 2 ) are the estimation of pitch and estimation error of pitch; (c 1 ) and (c 2 ) are the estimation of azimuth and estimation error of azimuth.Plotting the figure uses the triple Y-axis plotting approach because the results are in different orders of magnitude.

Figure 11 Figure 11 .
Figure 11 a 1 , a 2 , b 1 , b 2 , c 1 , c 2 illustrate the results of the north velocity, east velocity, and down velocity and the corresponding estimation errors of the MGWD system.The velocity estimated by PF, regardless of NEM or NNEM, apparently converges to zero velocity while the velocity error estimated

Figure 12 Figure 11 .
Figure 12 a1, a2, b1, b2, c1, c2 describe the estimation of latitude, longitude, and height of the MGWD system and the corresponding estimation errors.The latitude and longitude estimated by KF with LEM, PF with NEM, and PF with NNEM are almost close to the reference values 45.7776° for latitude, 126.6879° for longitude.The height estimated by PF with NNEM is 123.9996m (reference height: 124 m) while the error of height estimated by KF with LEM is about 5 m and the error of

Figure 12 a 1 ,
Figure 12 a 1 , a 2 , b 1 , b 2 , c 1 , c 2 describe the estimation of latitude, longitude, and height of the MGWD system and the corresponding estimation errors.The latitude and longitude estimated by KF with LEM, PF with NEM, and PF with NNEM are almost close to the reference values 45.7776 ˝for latitude, 126.6879 ˝for longitude.The height estimated by PF with NNEM is 123.9996m (reference height: 124 m) while the error of height estimated by KF with LEM is about 5 m and the error of height estimated by PF with NEM is approximately 3 m.The results demonstrate that the PF with NNEM can constrain the error growth of height effectively.Table 3 lists the results of Experiment 2.

Figure 12 .
Figure 12. Results of estimation of position and estimation error of position by KF with LEM, PF with NEM, and PF with NNEM under large-angle attitude error condition.(a1) and (a2) are the estimation of latitude and estimation error of latitude; (b1) and (b2) are the estimation of longitude and estimation error of longitude; (c1) and (c2) are the estimation of height and estimation error of height.

Figure 12 .
Figure 12. Results of estimation of position and estimation error of position by KF with LEM, PF with NEM, and PF with NNEM under large-angle attitude error condition.(a 1 ) and (a 2 ) are the estimation of latitude and estimation error of latitude; (b 1 ) and (b 2 ) are the estimation of longitude and estimation error of longitude; (c 1 ) and (c 2 ) are the estimation of height and estimation error of height.
Table 3 lists the results of Experiment 2.

Table 3 .
Navigation solutions and corresponding errors estimated by PF with NNEM, PF with NEM, and KF with LEM in Experiment 2.