Abstract
A circular error probability (CEP) metric in ballistic missile science is an experimental indicator of the accuracy of a missile system. There are a lot of error sources that cause a ballistic missile to deviate from its ideal trajectory, and that causes a deviation from required CEP. This work discusses the problems of dispersion of ballistic missiles due to inertial navigation system (INS) errors. INS deterministic errors are usually calibrated and compensated using some proper techniques. However, INS stochastic errors can be modeled and analyzed. In this study, a chosen missile is thoroughly analyzed using the six degrees-of-freedom missile flight trajectory simulator. A Monte Carlo simulation is used to generate a large number of flight trajectories to inspect the effect of INS stochastic noise on missile CEP. Moreover, a strategy for selecting an adequate sensor according to mission requirements and its corresponding sensor errors is introduced.
1 INTRODUCTION
Inertial navigation systems (INSs) are widely used in both military and commercial fields (Titterton et al., 2004). They are self-contained and reliable enough to satisfy users’ specifications in the most acceptable way possible. The most important and difficult ballistic missile issues (surface-to-surface) are guidance, navigation, and control problems. Guidance and control techniques are implemented to guide the missile through a reference path to reach the appropriate point of impact. The missile movement around the center of gravity (C/G) is achieved by the onboard autopilot, which helps in obtaining stability and transient restrictions. These guidance restrictions may be carried out by means of appropriate guidance and control techniques (Siouris, 2004).
One of these guidance techniques, mostly used with ballistic missiles, is inertial guidance. One of the most essential objectives of inertial guidance is to set the missile on a pre-calculated trajectory. The missile’s necessary attitude in various phases of flight is a function of the required range. There are many error sources that arise during missile guidance like INS errors, which affect the ability of inertial sensors to calculate the acceleration, position, and attitude of a missile.
Assisting the INS with external navigation sensors such as GPS may bound the propagating errors in INS algorithms. Nevertheless, this development is regulated by specific countries. For example, GPS receivers have certain limits on maximum height and velocity. Consequently, this work focuses on standalone inertial navigation systems. The precision of the navigation loop depends on the calibration, alignment, and algorithm of the adopted estimation techniques used with INS sensors (Zhang et al., 2012). A typical INS is comprised of inertial measurement units (IMUs) and the processing module that performs all the required navigation processing. The IMU in question is usually composed of two triads of gyroscopes (gyros) and accelerometers that measure angular rate and acceleration, respectively.
IMU errors are classified into systematic errors (deterministic) and random errors (stochastic). Calibration is performed for deterministic errors in accelerometers and gyros such as scale factor, non-orthogonality, bias, and misalignment. Calibration is done by going through different calibration techniques such as the six-position test, eight-position test, and angular rate test. The goal of IMU calibration is to estimate the intrinsic value of the IMU and compensate for errors in IMU measurements to enhance INS efficiency. The precision of the IMU intrinsic calibration will also have a direct impact on the effects of motion estimation (Hegazy et al., 2020a).
Stochastic noise analysis is used to estimate random errors such as angle random walk, bias instability, and rate random walk. Stochastic analysis is performed using several techniques such as the autocorrelation function, power spectral density, and Allan variance (AVR). The core of the analysis process is the calculation of the effect of specific noise sources on sensor output and sensor accuracy (Hegazy et al., 2020a).
A strap-down INS (SINS) typically contains an orthogonal three-axis group of gyros and accelerometers that provide data to the INS navigation computer. INS modules are directly installed on the vehicle chassis. Navigation is the process of calculating the physical body’s position, velocity, and acceleration relative to some reference frame (i.e., the inertial frame, Earth’s frame, navigation frame, and the wander azimuth navigation frame). The basic principle of SINS that is addressed in this work was calculated in the wander azimuth navigation frame (w-frame). The mathematical model of SINS was set up and the calculation steps in the w-frame were obtained (Zhang et al., 2012).
INS errors affect missile accuracy. A circular error probability (CEP) in ballistic strategic science is an experimental indicator of the accuracy of a given missile system. CEP can be defined as an indicator of the delivery accuracy of a weapon system, used as a factor in determining probable damage to a target (Moran Jr., 1966).
This work aims to enhance the precision of an INS by using appropriate and precise inertial sensor calibration methods to estimate both deterministic and stochastic errors. The inertial sensors adopted in this research are under the micro-electromechanical systems (MEMS) technology. To study the effect of INS errors on missile navigation performance, a complete missile six-degrees-of-freedom (6DOF) simulator is used. The INS is modeled according to the required noise parameters integrated with the missile model.
Missile performance is demonstrated by illustrating the difference between nominal trajectory in the onboard navigation computer and proposed algorithm outputs. The proposed algorithm contains an error model to subtract deterministic errors and an ideal gyro and accelerometer model to add stochastic noise to calculate the effect of stochastic noise on CEP using Monte Carlo simulations. The general model and subsequent algorithms were developed in the MATLAB environment. The results are demonstrated through simulation and real-flight trajectory data with an illustration of the effect of the IMU noise parameter on CEP. This workflow is illustrated in Figure 1.
CEP calculation according to selected IMU errors
The article layout is organized as follows: Section 1 focuses on the introduction; Section 2 regards INS sensor calibration, stochastic simulation, and analysis; Section 3 focuses on strap-down algorithms and Monte Carlo simulation; Section 4 describes the CEP algorithm; Section 5 regards experimental work and analysis; and Section 6 concludes the paper.
2 INS SENSOR CALIBRATION AND STOCHASTIC ANALYSIS
The main inertial sensor errors explored are predictable errors and unpredictable errors. Predictable (deterministic) errors are those that can be estimated and compensated for using a suitable calibration method, such as bias, scale factor, non-orthogonality errors, and gravity (g)-sensitive and non-g-sensitive drifts. Unpredictable (random) errors are correlated with sensor noise such as rate random walk, quantization noise, and bias instability, which can be analyzed and estimated using many techniques, among them the AVR technique which is adopted in this work.
2.1 IMU Calibration
An IMU is composed of three accelerometers and three gyroscopes. The following subsections discuss the methodology used to calibrate the used gyros and accelerometers for deterministic errors.
2.1.1 Gyroscope Calibration
A mathematical gyro error model of deterministic error can be obtained by:
1
where Wxo/p represents the gyro output of axis X; Kx represents the gyro scale factor; D0x represents constant drift terms; Wxi/p represents applied rate; Dyx and Dzx represent installation error coefficients of gyro axes Y and Z with axis X, respectively. Dx1, Dx2, and Dx3 represent drift coefficients of the X channel output of the gyroscope related to apparent acceleration in directions X, Y, Z, respectively; Wy and Wz denote the two components of angular rate respectively; and Ax, Ay, and Az represent applied input acceleration on the X, Y, and Z axes.
The gyro error model parameters were calibrated by making a relation between the sensor outputs with known inputs. The rate test alongside the multi-position static test is the most commonly utilized technique for calibration.
2.1.1.1 Scale Factor
Scale factor errors are analyzed by means of angular velocity calculations, in which the turntable includes a set of specific angular velocities for the under-test gyro axis. The scale factor is determined with the assumption that the gyroscope output has a linear relation to the input angular rate. Consequently, a constant scale factor is calculated as shown in Figure 2.
The scale factor’s trend of Axis Z with the angular input
However, this is not accurate for gyros since gyro damping non-linearly changes with varying input rates. This non-linearity problem is addressed in Guo and Zhong (2013) in which the authors propose a technique to evaluate scale factor taking into consideration the nonlinearity problem. The proposed technique calculates the scale factor in both directions (+ve scale factor and −ve scale factor) for each input rate, as shown in Figure 3. It appears that the scale factors have a linear relationship with the input, except for low rates (1°/s to 10°/s) due to the non-linearity.
The scale factor’s trend of Axis X with the input angular rates as positive and negative rates
Consequently, using standard processing techniques has the downside of obtaining different scale factors according to each input rate (+ve, –ve) which, consequently, provides unreliable scale factor estimation due to a wide range of interpolation and extrapolation. So, this approach uses the least-squares estimation as shown in Figure 4, depicting the relationship between the scale factors of Axis X with respect to the inverse of input angular velocity [−1 to – 30]−1 (° / s)−1 and [30 to 1]−1 (° / s)−1. This has drawbacks due to neglecting the scale factors at low rates. To overcome this downside, a calibration approach is proposed by this work’s author and presented in Hegazy et al. (2020a).
The least-squares curve for the regression of the scale factors
This innovative approach addresses the non-linearity problem of scale factor due to the shift in gyro damping with different rates applied in which a weighted linear regression is utilized to fit the output-input gyro response. Where the normalized error at each applied rate is calculated as shown in Figure 5.
The error behavior between the measured and actual rate
The weights are experimentally calculated as the normalized error inverse as follows:
2
where error(i) represents the scaled error rates by normalizing the error rates within a precise range [0, 1], and Hi represents the determined weight.
So, the scale factor calculated as:
3
where Kx,y,z is the single scale factor for the x, y, and z axes, and Ki is the scale factor for each rate.
2.1.1.2 Non-Orthogonality
Non-orthogonality errors occur when either of the sensor triad’s axes deviates from reciprocal orthogonality, usually occurring during the triad’s inertial sensor industrialization. There are many techniques to calculate the non-orthogonality error, such as the technique presented in Shin and El-Sheimy (2012), which proposes a non-orthogonality estimation algorithm from the static position test. This algorithm has a downside, however, in that it takes the Earth’s rotation as a single parameter, which is a weak signal. Moreover, it neglected the effect of centripetal acceleration from the lateral coupling factor (non-orthogonality).
In light of this problem, the work presented in Hegazy et al. (2020a) introduced an approach for the calculation of a pure non-orthogonal error by removing the irrelevant lateral coupling signal that is measured using a lateral accelerometer in position and rate tests. This irrelevant signal is directly proportional to the centripetal acceleration and is caused by the rates applied. The exact non-orthogonality error between y and x gyro is determined as follows:
4
where lyxT is the lateral coupling between Axis Y with respect to Axis X; Ayxg is the acceleration component in the y-axis, which identifies the gyroscope’s g-dependent drift; and Ky and K0y denote the scale factor of Axis Y’s gyroscope and accelerometer, respectively.
For the sake of verifying the proposal’s efficiency through a laboratory rate test, two error models are built from the classical approach presented in Guo and Zhong (2013) and in Shin and El-Sheimy (2012), for scale factor and non-orthogonality, respectively, and the proposed method. The turntable is utilized to apply different rates with the X-axis pointing up and the lateral Y-axis. The error of each error model with respect to the actual applied rate is depicted in Figure 6.
The error of the classical calibration method and proposed method with respect to the actual input rate
As shown in Figure 6, the proposed calibration method has a smaller error in both clockwise and counter-clockwise rates than the classical method. Moreover, it can be inferred that the difference between the two methods increases with high applied rates, which indicates the efficiency of the proposed method for high rate navigation systems.
2.1.2 Accelerometer Calibration
The multi-position static test is usually used for accelerometer calibration. Aggarwal et al. (2008) suggest that it is possible to construct an error calibration model by installing an IMU on a leveled turntable and alternately aligning each axis up and down for various angle rotations and measuring the scale factor and bias factor. However, this approach does not quantify the error of non-orthogonality. The modified six-position approach is then used to measure the IMU’s misalignment (non-orthogonality). Hegazy et al. (2020b) proposed a simple and more accurate calibration method that used an eight-position test to calculate all accelerometer error parameters without relying only on the direction of gravity acceleration. The position data is tabulated in Table 1.
Accelerometer Eight-Position Test
2.1.2.1 Scale Factor
While the input axis (IA) is up and down, two values are determined from the X accelerometer measurements using the previously explained sequence of positions for the accelerometer in the sensor assembly X direction, as follows:
Rotate around the output axis (OA) from Position 2 and Position 4:
5Rotate around the pendulum axis (PA) from Position 5 and Position 7:
6Then the scale factor is calculated from the average of these two values from Equation (5) and Equation (6):
7
2.1.2.2 Bias
For the accelerometer in the X direction, four values are calculated from the measurements of the X accelerometer using the previous explained series of positions while IA is up, east, down, and west as follows.
Rotate around the OA from Positions 3 and 1 and Positions 2 and 4:
8
9Rotate around the PA from Positions 6 and 8 and Positions 5 and 7:
10
11Hence, the bias is calculated from the average of these four values from Equation (8) and Equation (11):
12
2.1.2.3 Installation Error Coefficient
PA and OA misalignments are considered only because misalignments on the other IAs do not have a material impact on the accuracy of a single-axis accelerometer.
Rotate around the OA from Positions 1 and 3:
13Rotate around the PA from Positions 8 and 6:
14
Deterministic errors are compensated for using the error model in the INS algorithm after performing sensor calibration in the laboratory. However, stochastic noise, which affects INS output, is another cause of errors in INS sensors. As a result, accurate performance necessitates an investigation and analysis of these errors and their effects on the system, allowing an assessment of the output to choose the applicable sensors for any mission.
2.2 IMU Stochastic Simulation and Analysis
Gyroscope and accelerometer noise modeling is required for navigation simulation. To model such sensors, it is required to precisely study the noise characteristics such as quantization noise (Q), rate random walk (K), bias instability (B), and angle random walk (N; Hegazy et al., 2020c). IMU modeling is also used to assess the performance of error estimation algorithms such as the Allan variance (AVR) by calculating the previously modeled stochastic noise coefficients and, hence, the error percentage is able to be determined.
2.2.1 IMU Simulation
IMU simulation of stochastic noise is achieved using the MATLAB IMU Sensor System Object™. The IMU Sensor Object includes an idealized gyroscope and accelerometer object that is defined using certain parameters. First, using the noise-free default parameters for the gyro and accelerometer models, ideal sensor measurements are simulated. Second, the IMU’s parameters vary independently, since each parameter reflects a certain form of noise. The values of each parameter are determined with respect to the tactical grade sensor’s specification range.
The IMU simulation with different stochastic noise parameters was proposed by this work’s author and presented in Hegazy et al. (2020c). In brief, an example is presented here, such as how quantization noise (Q) has been introduced as a consequence of converting an analog signal to a digital signal. This is due to variations in the actual amplitudes of the points measured and the resolution of the analog-digital converter. The Q is generated as a white additive noise uniformly distributed within its bandwidth, and the simulation is achieved using the resolution parameter in gyro parameters given by the unit (rad/s)/LSB, which calculates the quantum step size. Figure 7 depicts a simulation in which the quantization error is measured as:
15
where WO/p denotes gyro output (quantization data) and Wi/p denotes angular velocity input.
Gyroscope input sinusoidal data and quantization data
Angle random walk (N) simulation is achieved using the noise density parameter, which is the amount of white noise in the sensor measurement; bias instability (B) is performed using the bias instability parameter representing the amount of pink or flicker noise in the sensor measurement; and the rate random walk (K) is simulated using the random walk parameter, which represents the amount of Brownian noise in the sensor measurement.
2.2.2 IMU Stochastic Analysis
The AVR technique is used to categorize various sources of noise found in measurements of inertial sensors. In this work, AVR is implemented and used in addition to the usage of AVR’s built-in MATLAB function. Validation of the AVR algorithm was performed by comparing the simulated input (added) noise to the estimated noise parameters and then calculating the error percentage.
The Allan variance is plotted in a log-log graph as the standard deviation versus the cluster time τ to estimate the IMU Crossbow 440 (Crossbow, 2010) noise parameters at different sampling frequencies after calibration and the removal of deterministic errors. Figure 8 shows the estimated noise parameters of Gyro X for four hours of data collected at a sampling frequency of 50 Hz. Gyro noise parameters at a sampling frequency of 50 Hz are tabulated in Table 2. The noise parameters of accelerometers are obtained and tabulated in Table 3.
The AVR analysis for gyro (IMU Crossbow 440) with sampling frequency of 50 Hz
Noise Parameter for Memsic-440 with Sampling Frequency of 50 Hz
Noise Parameter for Memsic-440 with Sampling Frequency of 50 Hz
3 MATHEMATICAL MODEL AND MATLAB SIMULATION OF STRAP-DOWN INERTIAL NAVIGATION SYSTEM
In a strap-down INS (SINS) system, the accelerometers and gyroscopes are attached to the vehicle body and don’t mechanically move. A mathematical approach is used to keep track of the IMU orientation and transfer the measurements from the body coordinate frame to the navigation coordinate frame to overcome the drawbacks faced by the gimbaled system and, most significantly, decrease the system’s scale, price, energy consumption, and complexity (Yin et al., 2013). The aim of the simulation of the trajectory was to produce data according to the designed real trajectory.
The mechanization of the wander azimuth is used to avoid all Earth navigation challenges (i.e., angular rotation rates and singularity issues; Jwo et al., 2014). The ENU-frame is selected as the frame for navigation. Figure 9 displays the whole SINS process principle in the wander azimuth frame mechanization.
Strap-down computation cycles
In this mechanization, the wander azimuth frame is defined as (xpypzp), in which the xp-axis points in the east direction; the yp-axis points in the north direction with an azimuth angle α; and the zp-axis points in the up direction. We define the Earth frame as (xeyeze) in which ye points to the North Pole, ze lies in the equatorial plane through the Greenwich meridian, and xe forms the right-hand coordinate frame. We define the navigation frame as (xn yn zn) in which xn points in the east direction, yn points in the north direction, and zn points in the up direction. Lastly, we define the body frame as (xb yb zb) in which xb points to right-wing, yb is the longitudinal axis, and zb points to the vertical (“up”) direction.
In order to define the orientation of the local-level frame in the wander azimuth system, a set of direction cosine matrices (DCMs) relating to Earth-fixed axes are utilized. It’s used to define the transformation matrix 
that is obtained by a series of rotations as shown in Figure 10.
Transformation cycle of reference frames for INS
If we define the platform frame as a wander azimuth frame or computational frame, it means that 
. Note that the platform frame (w-frame) is the same as the navigation frame (n-frame) when the wander azimuth angle equals zero.
The quaternion attitude is widely used in INS computation and is used to compute an equivalent DCM or transform the measured specific force vector into the chosen reference frame (Zhang et al., 2012).
First, the real trajectory (nominal trajectory) of the missile in the w-frame is set. Second, the ideal output of the IMU sensors is obtained from the nominal trajectory. The sensor simulation data can be acquired by adding noise (Crossbow IMU noise parameters) to the nominal data. Finally, we take the simulated sensor data with noise to extract an oblique simulated trajectory with stochastic noise parameters as shown in Figure 11, Figure 12, and Figure 13.
Accelerometer X output data from nominal trajectory versus data corrupted with noise
Gyroscope X output data from nominal trajectory versus data corrupted with noise
The error between nominal trajectory east and north distance and the trajectory corrupted with noise
Figure 11 and Figure 12 show the sensor data output of the nominal trajectory and the same data with added Crossbow noise parameters. The error between the nominal trajectory and the simulated trajectory corrupted with noise is shown in Figure 13. Calculation of the effect of each type of noise parameter on CEP requires a useful statistical tool for analyzing uncertain scenarios for noise parameters; such a method is called a Monte Carlo simulation (MCS).
4 MONTE CARLO SIMULATION
The MCS method is a common name for a wide variety of stochastic techniques. MCS is a very useful statistical tool for analyzing uncertain scenarios and provides a deep analysis of multiple scenarios. Information gathered from random samples is used to estimate the distributions and obtain statistical properties for different situations (Smits et al., 2018).
Monte Carlo simulation is a methodology used to research how a model responds to inputs that are randomly generated. Typically, it requires the following procedure:
Generate L inputs randomly (noise parameters).
For each of the L inputs, run a simulation. On a computerized model of the system being studied, simulations would be performed (according to the IMU model).
Accumulate and evaluate the results of the simulation. Popular measurements include mean value, distribution, and the minimum or maximum value of the output.
The sample distributions of the mean and the variance are close to normal distributions, as depicted in Figure 14 which shows the rate random walk sample distribution.
Normal Gaussian distribution of rate random walk
5 CIRCULAR ERROR PROBABILITY (CEP) ALGORITHM
The CEP circle is a statistical representation of the accuracy of a military strike; this index is proportional to the radius of a circle whose middle is at the target spot, and 50 percent of hits reach within this circle. The normal circle distribution is used to measure the CEP (Liu et al., 2018). Similar to statistical and probability analysis numbers, CEP fulfills:
16
where:
x denotes the horizontal range miss distance;
y denotes the vertical range miss distance;
R50 denotes CEP confidence levels; and
ρ denotes correlation coefficient.
The related bi-normal density distribution function is:
17
18
where σx, σy, μx, and μy represent the standard deviations and the means of horizontal and vertical range miss distance, and ρ is the correlation coefficient of X and Y. Thus, the CEP is determined by this integral equation.
Ballistic missile engineering has provided conclusive evidence that it is impossible for the horizontal miss distance to have any effect on the vertical miss distance and vice-versa. So, it is safe to assume that ρ = 0 (Didonato, 2007).
Many techniques can be used to resolve this equation and calculate the boundary of CEP, where the CEP is a circle considering the aiming point as its center and its diameter is the approximate value of the CEP. To solve this problem, a numerical method must be used.
The ideal method of approximating CEP using numerical integration of its bivariate normal density function is an infinite series expansion. The corresponding mathematical model for calculating the CEP is based on the normal modified distribution for the missile assuming the coordinate system is centered on the target with (0,0), and the deviation of long-range and cross-range is denoted as X and Y, respectively. The double integration for the joint probability density function is taken, and the boundary of this integration is the set of coordinates for a circle. This circle is taken as the CEP (Didonato, 2007). The equation for CEP determination is given by Equation (18) with substitute p= 0
19
A numerical method was used to solve this problem. We chose Simpson’s rule because it gives the minimum error during the evaluation process. Simpson’s rule is based on approximating curves with parabolas instead of line segments. The shaded area under the parabola is given by:
20
Applying this formula successively along a continuous curve Y = f (x) from x = a to x = b leads to an estimate of 
that is generally more accurate than the usual trapezoidal method.
By transforming Equation (19) to polar coordinates, with x = r cos θ and x = r sin θ yields:
21
and using elementary trigonometric identities:
22
where:
Using the fact that:
23
where IO(x) represents a modified Bessel function of the first kind with a zero-order consequently:
24
So:
25
6 EXPERIMENTAL WORK AND ANALYSIS
The experimental work structure is given in Figure 15. As shown in the figure, the gyroscope and accelerometer data are obtained from a nominal trajectory. First, the gyro and accelerometer data are corrupted individually and in combination with Crossbow stochastic noise parameters using an IMU modeling algorithm. For each noise parameter, a Monte Carlo simulation was used to produce 2,500 IMU data to produce 2,500 different trajectories (Liu et al., 2018). Finally, using the strap-down algorithm, the IMU data with its initial parameters are processed to obtain the system’s attitude, velocity, and position values. These position results yield 2,500 impact points, which are then used to calculate the CEP based on the deviation of each position in relation to the nominal trajectory position.
Simulation cycle
The CEP calculation is shown in Figure 16. As shown in the figure, the Crossbow IMU noise parameters caused the missile’s final impact point to be outside the required CEP (by 150 m in this case).
CEP with Crossbow noise parameters
So, a range for each noise type was selected in the scope of high tactical grade MEMS, simulated with gyro and accelerometer data to put hands on the permissible error ranges that lead to acceptable CEP values.
This is done for each noise individually and combined in the form of cases. Three such cases are mentioned as follows: Case 1 features the gyro angle random walk noise (N) parameter; Case 2 features the accelerometer’s N noise parameter; and Case 3 features the combined accelerometer and gyro N noise parameters. The results for these cases are illustrated in Figure 17, Figure 18, and Figure 19, respectively.
CEP with gyro angle random walk noise
CEP with accelerometer angle random walk noise
CEP N gyro and accelerometer
Consequently, the gyro and accelerometer noise parameters were also selected in the scope of high tactical grade MEMS and simulated with gyro and accelerometer data to put hands on the permissible error ranges that lead to acceptable CEP values. These parameters are listed in Table 4, and the CEP is shown in Figure 20.
Selection Noise Parameters
CEP for all IMU stochastic noise
7 CONCLUSION
The CEP circle is a statistical representation of the accuracy of a weapon system, so the paper has presented a method for the improvement of CEP through precise identification, calibration, and compensation of sensor deterministic error coefficients.
Precise modeling and simulation of different noise parameters were discussed, and the Allan variance algorithm was performed and validated. Moreover, the MCS method was introduced as a statistical tool for analyzing uncertain scenarios and providing a deep analysis of different scenarios for stochastic noise parameters. Also, a navigation and CEP algorithm was built for measuring the CEP, which indicates the ballistic missile accuracy, and estimates the effect of each stochastic noise parameter on CEP to choose sensors that are adapted to various applications.
HOW TO CITE THIS ARTICLE
Hegazy, S., Kamel, A., Arafa, I., & Elhalwagy, Y. (2022). INS stochastic noise impact on circular error probability of ballistic missiles. NAVIGATION, 69(2). https://doi.org/10.33012/navi.523
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.