## 3 Methodology
### MDDA Augmented GNSS/INS Integration
#### 3.6.2 INS Mechanization  
Local level frame (LLF) mechanization equations have been requested in many applications because the user may easily understand the navigation solution provided by the LLF navigation equations when they are on or near the surface of the Earth [49]. In the ECEF frame, the position vector \(p^{l}\) is defined as shown in Equation (21):  
\[p^{l}= \begin{bmatrix}\varphi\\ \lambda\\ alt\end{bmatrix}  \tag{21}\]  
where \(\lambda\) is the longitude \(long\), \(\varphi\) is latitude \(lat\), and \(alt\) is altitude. The velocity components influence the temporal rate of variation of the position components [3]:  
\[ \begin{bmatrix}\dot{\varphi}\\ \dot{\lambda}\\ alt\end{bmatrix} = \begin{bmatrix}\frac{1}{R_{m}+alt}&0&0\\ 0&\frac{1}{\cos\varphi(R_{m}+alt)}&0\\ 0&0&-1\end{bmatrix}\begin{bmatrix}vel_{n}\\vel_{e}\\vel_{d}\end{bmatrix} ^{l} \tag{22}\]  
\[\dot{p}^{l}=A^{-1}vel^{l} \tag{23}\]  
where \(\lambda\), \(\dot{\varphi}\), and \(alt\) are temporal rate of variations of long, lat, and alt. \(vel_{n}\), \(vel_{e}\) and \(vel_{d}\) are the velocities in the north, east, and down (NED) directions, respectively. \(R_{n}\) and \(R_{m}\) are normal and meridian ellipsoid radii. \(A^{-1}\) is the transformation matrix of the velocity vector between rectangular and curvilinear coordinates. The time derivative of LLF velocity component \(vel^{l}\) is shown in Equation (24). \[vel^{l}=C_{b}^{l}ax^{b}-(2\Omega_{ie}^{l}+\Omega_{el}^{l})vel^{l}+g^{l} \tag{24}\]  
where the \(C_{b}^{l}\) matrix transforms the b-frame to the l-frame, \(ax^{b}\) is the accelerometer's measurement in the b-frame, \(g^{l}\) is the vector of gravity in the l-frame, and \(\Omega_{ie}^{l}\), \(\Omega_{el}^{l}\) are the skew-symmetric matrices related to \(\omega_{ie}^{l}\), \(\omega_{el}^{l}\), which are Earth rotation rate and transportation rate [3]. The orientation (attitude) is calculated by solving the transformation matrix's time derivative equation \(\dot{C}_{b}^{l}\) as follows:  
\[\dot{C}_{b}^{l}=C_{b}^{l}[\Omega_{ib}^{b}-C_{l}^{b}\;(\Omega_{ie}^{l}+\Omega_{ el}^{l})C_{b}^{l}] \tag{25}\]  
where \(\Omega_{ib}^{b}\) is the skew-symmetric matrix of the gyro-measurement vector. The quaternion was used for the following reasons [3; 49; 50]:  
1. The singularity issue that can occur when using Euler angles is avoided by the quaternion solution. \begin{table}
\begin{tabular}{l c c} \hline \hline  & **Gyroscope** & **Accelerometer** \\ \hline Max measurement range & \(\pm 2000\) deg/s & \(\pm 8\) g \\ \hline Bias error & \(\pm 0.4\) deg/s & \(\pm 0.6\) mg \\ \hline Output noise & \(0.007\) deg/s/\(\sqrt{\text{Hz}}\) & \(180\) g/\(\sqrt{\text{Hz}}\) \\ \hline Resolution & \(4.8\times 10^{-4}\) deg/s & \(\pm 0.975\) mg \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Specifications of MEMS inertial sensors for the used smartphone [48]. It is rather easy to execute the quaternion computation. 3. Instead of six differential equations if directions cosines are used, just four equations are numerically solved. Quaternion \(q\) has four parameters \((q_{1},q_{2},q_{3},\) and \(q_{4})\) and is represented as follows: \[q=\left[\cos(\frac{\theta}{2}),\ (\frac{\theta_{x}}{\theta})\sin(\frac{\theta}{2}), \ (\frac{\theta_{y}}{\theta})\sin(\frac{\theta}{2}),\ (\frac{\theta_{z}}{\theta})\sin(\frac{ \theta}{2})\right]^{T}\] (26) where \(\theta\) is the rotational angle of the quaternion and equals \(\sqrt{(\theta_{x}^{2}+\theta_{y}^{2}+\theta_{z}^{2})}\) and \(\frac{\theta_{x}}{\theta},\ \frac{\theta_{y}}{\theta},\ \frac{\theta_{z}}{\theta}\) is the directional cosine of the rotation axis. Additionally, quaternion obeys the unity restriction where \(1={q_{1}}^{2}+{q_{2}}^{2}+{q_{3}}^{2}+{q_{4}}^{2}\) Hence, \[\vartheta=\frac{1}{2}[\tilde{\Omega}(\omega)]q\] (27) \[\tilde{\Omega}(\omega)= \begin{bmatrix}0&\omega_{z}&-\omega_{y}&\omega_{x}\\ -\omega_{z}&0&\omega_{x}&\omega_{y}\\ \omega_{y}&-\omega_{x}&0&\omega_{z}\\ -\omega_{x}&-\omega_{y}&-\omega_{z}&0\end{bmatrix}