Indoor Positioning Integrating PDR/Geomagnetic Positioning Based on the Genetic-Particle Filter

Abstract

This paper proposes a fusion indoor positioning method that integrates the pedestrian dead-reckoning (PDR) and geomagnetic positioning by using the genetic-particle filter (GPF) algorithm. In the PDR module, the Mahony complementary filter (MCF) algorithm is adopted to estimate the heading angles. To improve geomagnetic positioning accuracy and geomagnetic fingerprint specificity, the geomagnetic multi-features positioning algorithm is devised and five geomagnetic features are extracted as the single-point fingerprint by transforming the magnetic field data into the geographic coordinate system (GCS). Then, an optimization mechanism is designed by using gene mutation and the method of reconstructing a particle set to ameliorate the particle degradation problem in the GPF algorithm, which is used for fusion positioning. Several experiments are conducted to evaluate the performance of the proposed methods. The experiment results show that the average positioning error of the proposed method is 1.72 m and the root mean square error (RMSE) is 1.89 m. The positioning precision and stability are improved compared with the PDR method, geomagnetic positioning, and the fusion-positioning method based on the classic particle filter (PF).

1. Introduction

With the development of modern society, people’s demands for location services is becoming increasingly urgent [1]. The Global Navigation Satellite System (GNSS) [2] can provide users with a high-precision positioning service in the outdoor environment. However, in the indoor environment, because the GNSS signal is blocked by humans, buildings and many other factors, it becomes so weak that the positioning accuracy will be greatly reduced. In order to solve indoor positioning problems, many kinds of method have been studied, such as Wi-Fi positioning [3,4], Bluetooth positioning [5,6], ultra-wideband (UWB) [7,8], radio frequency identification (RFID) [9], infrared [10], ZigBee [11], etc. Every positioning method is capable of providing an indoor location service for people, especially the UWB technology which can reach sub-meter accuracy in the indoor open areas [8]. But these approaches require additional signal transmitting devices, which means that people cannot obtain their position if there is no such equipment in the environment. Moreover, these devices are usually expensive which will make the cost very high if we apply these technologies to a large-scale environment.

Pedestrian dead-reckoning (PDR) [12,13] technology can use an accelerometer, gyroscope, and magnetometer to obtain continuous positions instead of installing the signal transmitting stations. But the positioning error will accumulate quickly because of the initial position error, heading angle estimation error and other measurement errors. Therefore, it is necessary to integrate PDR with other positioning methods to restrict error accumulation.

In 2000, Suksakulchai et al. [14] conducted an experimental analysis in the corridor environment and proposed that the indoor magnetic field can be used for indoor positioning. Then, experiments in [15] studied in detail the stability and distribution of an indoor magnetic field, pointing out that it is stable and reliable. All of these related research works provide a solid theoretical basis for the application of the magnetic field for indoor positioning. In recent years, researchers have found that utilizing indoor magnetic field can achieve high positioning precision at almost no investment that are safe and infrastructure-free compared with other technologies [16]. These characteristics make indoor geomagnetic positioning has great research value and also gradually become the research hotspot.

Today there are many kinds of geomagnetic positioning methods, which can be divided into two main types: sequence-based magnetic matching positioning (SBMP) and single-point-based magnetic matching positioning (SPMP). The dynamic time-wrapping algorithm (DTW) [17] is a representative positioning algorithm of SBMP. In [18,19], experimenters applied the DTW algorithm to magnetic matching by comparing the similarity of different trajectories and employed this magnetic positioning method on smartphones. The results show that the accuracy of 80% test positions is approximately 2 m. Generally speaking, the SBMP method has high positioning precision, but it relies on stable magnetic data sequences. However, the magnetic data sequence will generate large fluctuations due to the device heterogeneity. Even though some methods have been studied to solve this problem, they cannot be implemented in real-time positioning [20,21,22]. What is worse is that SBMP is restricted by the trajectory and walking speed [16,23]. In other words, if participants want to get their position, they should walk on the trajectory where they collect magnetic field data before, which means that this method will perform poorly in the open area.

By contrast with SBMP, the SPMP method uses the geomagnetic features data at the signal points rather than the magnetic sequence data and there is also no limitation to the speed and trajectory of pedestrian walking. Therefore, it is more flexible and easier to implement compared with SBMP method. This method usually does not have high positioning precision if there are few geomagnetic features. But the particle filter algorithm makes up for this shortcoming and has high positioning accuracy, which is also a typical single-point-based fingerprint positioning algorithm [24,25]. In related research, Huang et al. [26] used the Hausdorff distance and Pearson correlation coefficient [27] to control the initial position error and accelerate the convergence speed of the filter. In [28], an adaptive particle filter algorithm was proposed in geomagnetic positioning. The number of particles in the experiment can adjust according to the different features and the average positioning accuracy was 0.51 m. Shu et al. [29] combined the Wi-Fi signal with magnetic field data and designed a two-pass bidirectional particle-filtering process to improve precision, which achieved a mean location accuracy of 0.9 m in the office environment. Although geomagnetic positioning technology based on the particle filter can achieve high accuracy, it suffers from the particle degradation problem. Increasing particle numbers is one of the approaches to solve this problem, for example, the number of particles in [28] even reaches 3000. A large number of particles will result in a large computational load and limit its application on smart terminals. Therefore, how to better solve the particle degradation problem is the key to improving the performance of the particle filter.

In this paper, we choose the SPMP method for geomagnetic positioning because it is more flexible than the SBMP method. Five magnetic features data are extracted to ensure the reliability of geomagnetic fingerprint and the accuracy of the geomagnetic multi-features positioning algorithm. Moreover, the Mahony complementary filter (MCF) algorithm is applied to calculate the heading angles in PDR technology. In order to optimize the particle degradation problem, we use the gene mutation method in the genetic algorithm [30] and the method of reconstructing the particle set by combining the geomagnetic positioning results to update the particle set. Lastly, we integrate the PDR and the geomagnetic positioning by using the genetic-particle filter (GPF) algorithm and carry out several experiments to evaluate its performance. The contributions can be summarized as follows:

• We transform the geomagnetic data into the geographic coordinate system (GCS) and extract five geomagnetic features data. These features data can improve the specificity of geomagnetic fingerprint and the accuracy of the devised geomagnetic multi-features positioning algorithm.
• We design the optimization mechanism for the particle degradation problem, which utilizes the genetic mutation method to increase the diversity of the resampled particles and uses the particle set reconstruction method by combining the geomagnetic positioning results to improve the particles’ reliability. The proposed mechanism better ameliorates the particle degradation problem;
• We propose the fusion-positioning method that utilizes the PDR positions as the system state and uses the geomagnetic positions as the system observation based on the genetic-particle filter. The mean positioning error is 1.72 m and the root mean square error is 1.89 m.

The remainder of the paper is organized as follows: Section 2 gives a detailed description of each module, such as PDR technology; Section 3 evaluates the performance of the proposed methods; Section 4 summarizes the conclusions and discussion.

2. Materials and Methods 2.1. Pedestrian Dead-Reckoning (PDR) Module

PDR technology mainly includes three parts: heading angle estimation, step detection, and step length estimation. As illustrated in Figure 1, if the initial position is known, PDR technology can utilize the measured heading angles and step length to reckon positions by using Equation (1). However, due to the factors like the accuracy of heading estimation and the error of the initial position, system error will accumulate rapidly with time, which directly causes the navigation trajectory to deviate seriously from the actual path. Therefore, other positioning methods are usually needed to weaken the influence of error accumulation on the positioning results.

{_Nk+1=_Nk+_lk•cos_ψkEk+1=_Ek+_lk•sin_ψk

where_Nkand_Ekrepresent the coordinate values of the north and east at timek, respectively;_lkand_ψkrepresent the step length and heading angle at timek, whose calculation methods will be described in the following sections.

2.1.1. Heading Angle Estimation

Estimating heading angle is important for the PDR method to calculate position coordinates. In general, the heading angle can be measured by using the gyroscope. But there are many integration calculations that will make the measured heading angle inaccurate. In our research, the Mahony complementary filter [31] is implemented on the smartphone and used to compute the heading angle. It should be clear that the MCF algorithm mainly utilizes the gyroscope data to calculate the attitude of carrier, but the accelerometer and magnetometer will correct the gyroscope data so that the estimation error will not accumulate.

Generally, the whole producer of the MCF algorithm contains error compensation calculation and gyroscope data correction. The gyroscope data will be expressed as the quaternion vector, which represents current attitude prediction. Assuming that the quaternion vector is Q = [q₀ q₁ q₂ q₃]^T, it is defined as follows:

{_q0=cos(β/2)_q1=_ωx′sin(β/2)_q2=_ωy′sin(β/2)_q3=_ωz′sin(β/2)

whereβis the rotation angle,_ωx′,_ωy′and_ωz′are the normalized instantaneous gyroscope data. Then, in the error compensation calculation stage, two error items will be calculated by using the accelerometer and magnetometer data. If the gyroscope error correction ise=^[_exeyez]T, it can be defined as:

e=_ea+_em

where_eaand_emare the error correction items calculated by using the accelerometer and magnetometer data, respectively. They can be expressed as_ea=^{[_eaxeayeaz]T}and_em=^{[_emxemyemz]T}. Their calculation methods are as follows:

{_ea=(_Cnb•_ga)×a_em=(_Cnb•_bm)×m

where_gais the gravity vector in the geographic coordinate system,_ga=^[00g]T,_bmis the geomagnetic vector when x-axis of device points to the north,_bm=^[_bmx0_bmz]T;aandmare the normalized measured accelerometer and magnetometer data;_Cnb is the rotation matrix from the geographic coordinate system to the carrier coordinate system and it will be described in Section 2.2; “×” represents the vector cross product, respectively. After getting the error correction, the corrected gyroscope data can be calculated as:

ω=_ωg+_KPe+_KI∫e

where_ωg=^{[_ωgxωgyωgz]T}andω=^{[_ωxωyωz]T}are the normalized gyroscope raw data and the corrected gyroscope data, respectively;_KPand_KIare the error control items, which are 2.0 and 0.001 based on our test, respectively. After getting the corrected gyroscope data, we should substituteω into the quaternion differential equation and solve this equation by using the first-order Runge–Kutta method [32] to obtain the quaternion update equation as Equation (6). The quaternion values can be calculated by this equation:

{_q0(t+T)=_q0(t)−T2[_{ωx q1}(t)+_{ωy q2}(t)+_{ωz q3}(t)]_q1(t+T)=_q1(t)+T2[_{ωx q0}(t)+_{ωz q2}(t)−_{ωy q3}(t)]_q2(t+T)=_q2(t)+T2[_{ωy q0}(t)−_{ωz q1}(t)+_{ωx q3}(t)]_q3(t+T)=_q3(t)+T2[_{ωz q0}(t)+_{ωy q1}(t)−_{ωz q2}(t)]

where_qi(t)and_qi(t+T)are the quaternion values at timetand(t+T),i=1,2,3,⋯,T is the sampling time interval. Then, normalizing the updated quaternion values and substituting them into Equation (7), the heading angle can be obtained. In our research, the first-order low-pass filter [33] and the first-order high-pass filter [34] are applied to process the observation noise of the acceleration raw data and the gyroscope raw data, respectively. Because there are serval heading angles in one second, the median filter [35] is used to reduce the measurement error and obtain more accurate results.

{ψ^=tg−12(_{q1 q2}+_{q0 q3})1−2(_q22+_q32),ψ∈(0,2π)θ=−^sin−12(_{q1 q3}−_{q0 q2}),θ∈(−π2,π2)γ^=tg−12(_{q2 q3}+_{q0 q1})1−2(_q12+_q22),γ∈(−π,π)

whereψ,θ, andγare the heading angle, pitch, and roll, respectively,ψis also called yaw.

2.1.2. Step Detection

There are many step detection methods in PDR technology, such as peak detection [36], zero velocity update [37], etc. As shown in Figure 2a, the acceleration changes regularly during the walking process and the acceleration cure has peaks and valleys. The peak detection method realizes the step detection by detecting the peaks and valleys. This method is easy and has a low calculation burden. However, the measurement error will affect it when recognizing steps. In our research, we employ the finite-state machine (FSM) algorithm to the step detection process because this algorithm is easy to be implemented and more resistant to error interference than the peak detection method.

In fact, FSM is a special peak detection algorithm. It divides the walking process into five states, which can be recorded as_S0~ _S4._S0represents the steady state and the acceleration at this time is approximately equal to 9.8 m/s²;_S1is the taking step state and the acceleration gradually increases;_S2is the finding peaks state after entering the state_S1;_S3is the peak threshold judgment state and filter out the false peaks;_S4 represents the falling state and the acceleration gradually decreases. As shown in Figure 3, the current walking state is transferred in the above five states and one step will be recognized if the five states all satisfy the threshold conditions. When a large measurement noise appears during the algorithm execution, it will return to the_S0state and re-detect the walking state. Therefore, FSM can effectively eliminate the interference of false peaks and the measurement errors in the step detection process.

Table 1 gives all the variable names of FSM. It should be noted that_Pupand_Pdownwill be reset as zero at_S0state or there exist large measurement errors. They can be calculated as follows:

{_Pup=_Pup+1,_ai−_ai−1>acc_diff_Pdown=_Pdown+1,_ai−_ai−1<acc_diff

where_aiand_ai−1are the acceleration values at timeiandi−1, respectively.

2.1.3. Step Length Estimation

When people are walking, the step-length can be estimated by using the acceleration data. The existing step length models include linear model [38], constant models [39], nonlinear models [40], etc. Because of the difference between people’s height and stride, the step length is also different. Therefore, the application of the constant model and linear model has certain limitations. We use the Weinberg model [41] that is a non-linear model to estimate step length. During the one-step process, if the maximum and minimum of the acceleration are A_max and A_min, step length can be calculated as:

StepLength=K•_Amax−_Amin4

where K is a constant value, which can be obtained by training a large number of acceleration data.

2.2. Geomagnetic Positioning Module

2.2.1. Geomagnetic Multi-Features Data Extraction

It will be better to describe the distribution of the magnetic field at a single point by utilizing more geomagnetic features data. In addition to the three-axis components(_BxByBz), the horizontal intensity_Bhand the total intensityBrepresent the magnitude of the magnetic field in the horizontal plane and the 3-dimensional space. If these two components are utilized for geomagnetic matching positioning, they can provide additional distance information for the matching process. Therefore, we extracted these two features and their calculation methods are as follows:

{_Bh=_Bx2+_By2B=_Bx2+_By2+_Bz2

However, the three-axis components of the magnetic field measured by smartphones are in the carrier coordinate system (CCS), whose x, y and z axes point to the right, front and top of the device as shown in Figure 4a. These three components vary with the attitudes of smartphones and it is meaningless to directly use the untransformed three-axis magnetic data. Therefore, it is necessary for us to transform the three-axis geomagnetic data to the geographic coordinate system (GCS) because they are stable in the GCS.

The GCS and the CCS can be rotated into each other through attitude angles, which contain yaw, pitch, and roll. Their calculation methods have been described particularly in Section 2.1.1. The rotation matrix_Cnbbetween the two systems can be defined as:

_Cnb=_{Cn3 Cn2 Cn1}=[cosγ0−sinγ010sinγ0cosγ][1000cosθsinθ0−sinθcosθ][cosψ−sinψ0sinψcosψ0001]

where_Cn1,_Cn2and_Cn3 are the corresponding matrices for yaw, pitch, and roll. Yaw, pitch, and roll respectively represent rotation around the z, x, and y axes as shown in Figure 4. The rotation matrix_Cbnfrom the CCS to the GCS can be calculated as:

_Cbn=^(_Cnb)−1=^(_Cnb)T

where^(_Cnb)−1and^(_Cnb)Tare the inverse matrix and transposed matrix of_Cnb. Supposing that the three-axis components vector of magnetic data in the CCS is_BM=^[_BxByBz]Tand_BS=^{[_Bx ′_By ′_Bz ′]T}is the magnetic data vector in the GCS, the relationship between them is defined as:

_BS=^(_Cnb)T _BM=[cosγcosψ+sinγsinψsinθsinψcosθsinγcosψ−cosγsinψsinθ−cosγsinψ+sinγcosψsinθcosψcosθ−sinγsinψ−cosγcosψsinθ−sinγcosθsinθcosγcosθ][_BxByBz]

We can obtain three-axis stable geomagnetic data by using Equation (13). As shown in Figure 5, the transformed data are stable and the x-axis components data are close to zero, which are consistent with theoretical results. However, we found that the x-axis components data of some experimental areas were not close to zero in our research. Therefore, we retain the x-axis component in this paper. Eventually, we can extract five geomagnetic features data as the fingerprint data by using the above method.

2.2.2. Geomagnetic Multi-Features Positioning Algorithm

The geomagnetic multi-features positioning algorithm includes the offline database construction stage and online positioning stage. During the database construction stage, the experimenter uses mobile phones to collect the geomagnetic data on each selected point for 15 s, then extracts five geomagnetic features data. The fingerprint of a single point can be expressed as(_Bx,_By,_Bz,_Bh,B,x,y) by combining the actual geographic coordinates. For the points without data collection, a linear interpolation [42] is used to obtain the geomagnetic data of these points.

In the online positioning stage, the real-time collected geomagnetic data will be firstly transformed into the geographic coordinate system, and five geomagnetic features data will be extracted. Then the k-nearest neighbor (KNN) method [43] is employed to obtain the positioning results by traversing the geomagnetic database. If one piece of the real-time collected data is(_Bx ^′,_By ^′,_Bz ^′,_Bh ^′,^B′), the distance between the measured and fingerprint data can be calculated as:

d=^{(_Bx−_Bx ′)2}+^{(_By−_By ′)2}+^{(_Bz−_Bz ′)2}+^{(_Bh−_Bh ′)2}+^(B−B′)2

If there arenfingerprint data in the database,ndistances can be obtained by traversing the geomagnetic database. We should find thekminimum distances, in the meantime, the coordinate points corresponding to thesekminimum distances should be recorded. If the current position vector is_Xm=^[_xmxymy]T, it can be calculated as:

_Xm=^[_xmxymy]T=^{[1k•∑i=1k_xi1k•∑i=1k_yi]T}

where_xiand_yiare the coordinate values of theknearest points, i = 1,2, ⋯ k.

In one second,mpieces of geomagnetic data can be collected. All of the collected geomagnetic data will be involved in the positioning calculation. Thus, there will bem positioning results in one second. These results are essential for fusion positioning in Section 2.3. In order to better understand the geomagnetic multi-features positioning algorithm, Algorithm 1 gives the detailed process of it.

Algorithm 1 Geomagnetic multi-features positioning algorithm
1:	Input: magnetic database MD, Coordinate data CD, measured magnetic dataBmeasure, positive integerk
2:	Bmeasure←normalize(Bmeasure)
3:	fori← 1 to length(Bmeasure)do
4:	forj← 1 to length(MD) do
5:	calculate the distance dj
6:	add toDist,Dist← dj
7:	end
8:	findkminimum values of theDist,mindist←min(Dist)
9:	index←find(Dist==mindist)
10:	current_coordinate←sum(CD(index))/k
11:	addcurrent_coordinatetopos_result,pos_result←current_coordinate
12:	end
13:	Output:pos_result

2.3. Fusion Positioning Module

2.3.1. Classical Particle Filter

The particle filter (PF) is an algorithm based on the Monte Carlo method [44] and the Bayesian approach [45], which has good performance in solving non-linear problems. Generally, a dynamic system can be assumed that it obeys the Markov process of order one, it can be modeled as follows:

{_Xk=_fk(_Xk−1)+_δk−1Zk=_hk(_Xk)+_σk

whereX(k)andZ(k)are the system state and observation at timek,_δk−1and_σkare the process noise and measurement noise at timek−1andk, respectively._fkreflects the relationship between the current and previous state._hkexpresses the relationship between states and observations. The posterior probability density function (pdf)p(_Xk|_Z1:k)is usually used to calculate the optimal estimation and predict the appearance probability ofX(k), which is called as the prediction in PF.p(_Xk|_Z1:k)relies on the previous posterior pdfp(_Xk|_Z1:k−1)based on the Bayesian theory, which can be calculated by using Equation (17) if there are observations available:

P(_Xk|_Z1:k−1)=∫p(_Xk,_Xk−1|_Z1:k−1)d_Xk−1=∫p(_Xk|_Xk−1)p(_Xk−1|_Z1:k−1)d_Xk−1

wherep(_Xk|_Xk−1)is the transition probability distribution defined byX(k),p(_Xk−1|_Z1:k−1)is the previous posterior, respectively. Assuming the previous posterior is known, the current posterior is defined as follows:

P(_Xk|_Z1:k)=p(_Zk|_Xk,_Z1:k−1)p(_Xk|_Z1:k−1)p(_Zk|_Z1:k−1)=p(_Zk|_Xk)p(_Xk|_Z1:k−1)p(_Zk|_Z1:k−1)

wherep(_Zk|_Z1:k−1)is the normalizing constant, it is calculated by using Equation (19):

p(_Zk|_Z1:k−1)=∫p(_Zk|_Xk)p(_Xk|_Z1:k−1)d_Xk

We can use Equation (18) to update the current posterior probability and this is an important step in PF. But it is difficult for us to solve Equations (17) and (19) due to the complex integral operation. The Monte Carlo method provides an effective approach for solving this problem. This method uses a sequence of weighted samples to represent the posterior pdf and the mean value of these samples is substituted for the integration. If the sample set hasNsamples,p(_Xk|_Z1:k)can be calculated by:

p(_Xk|_Z1:k)≈1N_∑i=1Nωkiδ(_Xk−_Xki)

where_ωkiand_Xkiare the weight and the state of i-th sample,δ(⋅)is the Dirac delta function. The update formula of_ωkiis defined as follows:

_ωki∝_ωk−1ip(_Zk|_Xki)p(_Xki|_Xk−1i)q(_Xki|_Xk−1i,_Zk)

whereq(⋅)is the importance density. WhenN→∞ , the result of Equation (20) approaches the true posterior density [45]. PF will resample a different set at each timestampkbased on the updated posterior probability, which can also be called weight. Resampling will copy large-weight samples and abandon the small-weight ones to obtain the optimal estimation. But this will decrease the diversity of samples if the system runs for a long time. This phenomenon is particle degradation which can affect the performance of PF and waste computing resources. Therefore, ensuring the diversity of samples and optimizing particle degradation are the key factors in determining the performance of PF. We concentrate on these problems and design the optimization mechanism for particle degradation in our research work.

2.3.2. Gene Mutation

The gene mutation method of the genetic algorithm (GA) [30] is essential to genetic algorithm optimization. GA follows the “survival of the fittest” principle and it uses the genetic operators of natural genetics to cross and mutate the population to create new generations. The information of each individual is expressed as the binary code, which is called the DNA. The whole process of gene mutation is as shown in Figure 6.

Assuming that the information of an individual isIand it can be described in the decimal system. In the mutation process,Iwill be transformed into the binary code vectorM, which can be described as:

M=_{[bibi−1⋯b0]2}

where_biis the binary code, its value is 0 or 1. And then, the binary code will change randomly to imitate the DNA mutation and individuals with different characteristics can be obtained by decompiling the mutated binary codes using Equation (23):

I=_{(∑j=0ibj•^2j)10}

2.3.3. Optimization Mechanism for Classical Particle Filter

The optimization mechanism aims at improving the diversity and reliability of the resampled particles to ameliorate the particle degradation. Therefore, we utilize the gene mutation method to increase the diversity and use the reconstruction particle set method by combining the geomagnetic positioning results to improve the reliability of particles. We can assume that the characteristics of each resampled particle are its coordinate values. Firstly, we should utilize the gene mutation method to transform the coordinates into binary codes. Since the coordinate values are decimal, the problem becomes the decimal-to-binary conversion question [46]. As shown in Algorithm 2, it gives the detailed transformation process ofxvalue. It should be clear that we should outputMin the inverse order to obtain the binary codes.

Algorithm 2 decimal-to-binary ofxvalue
1:	Input: coordinate valuex, empty setM, Boolean variableisRun
2:	get the integer part ofx,x_in←floor(x);isRun← true
3:	while (isRun) do
4:	calculate the remainderr,r←x%2
5:	calculate the quotientq,q←(x−r)/2
6:	ifq==0do
7:	calculate the remainderr,r←x%2
8:	addrtoM,M←r;isRun←false
9:	else do
10:	addrtoM,M←r
11:	end
12:	x_in←q
13:	end
14:	binarycode← reverse (M)
15:	Output:binarycode

The binary codes ofy value can be obtained by using Algorithm 2. As shown in Figure 7, if we change the binary codes to imitate the gene mutation, we can get the particle swarms with different characteristics. Therefore, the resampled particles become various and the diversity is improved.

Apart from increasing the diversity, we should also assure the reliability of particles so that PF abandon fewer particles in the resampling step. In Section 2.2.2, we have pointed out that there are many geomagnetic positioning results rather than one result. In these results, we should selectsresults that are close to the PDR position. If the PDR positioning results areX(k)={_xk,_yk}and the geomagnetic positioning results areZ(k)={_Z1,_Z2,⋯,_Zm},_Z1=(_x1,_y1),⋯, the Euclidean distance betweenZ(k)andX(k)can be computed by using Equation (24):

_Di=^(_xi−_xk)2+^(_yi−_yk)2

Then, based on the system stateX(k) , we can also use the KNN method in Section 2.2.2 to findsgeomagnetic results that are close to the PDR position to replace thes particles. After this optimization, each generation particles are composed of the previous generation “excellent” particles after mutation and the geomagnetic positioning results of the current generation. These two methods assure the diversity and reliability of the particle set and avoid the “excellent” particles reduction of each generation. The particles of each generation seem to be in evolution rather than degradation in classical PF. Figure 8 shows the detailed process of the genetic-particle filter algorithm in this paper, and the size of different balls represents different characteristics of particles.

2.3.4. Fusion-Positioning Algorithm based on Genetic-Particle Filter

In this section, we will describe the detailed process of the fusion positioning algorithm based on the G-Particle Filter. We use PDR as the system stateX(k)and the geomagnetic positioning resultsZ(k)={_Z1,_Z2,⋯,_Zm},_Z1=(_x1,_y1),⋯, as the system observation; the PDR position will be corrected in real-time by the fusion positioning results. The general process is as follows:

Step 1, predict the system stateX(k)at timekand reconstruct the particle set based onX(k). We assume the particle set isP(k)={_P1,_P2,⋯,_Pn}and the weight of each particle is1/N. Then, we should select thesobservations with the smallest distance to replace the lastsparticles to reconstruct the particle set. The new particle set isP(k)={_P1,_P2,⋯_Pn−s ^′,⋯,_Pn ^′}.

Step 2, update and normalize the weights of particles. After getting the reconstructed particle set, each particle’s weight should be updated by using the Gaussian function:

_ωi(k)=1_σr2πexp[−^{‖_Pi−X(k)‖2}2_σr ²]

where_σris the observation noise. The normalized weight of each particle can be calculated by:

_^ω′i(k)=_ωi(k)∑i=1n_ωi(k)

where_ωi(k)and _^ω′i(k)are the un-normalized and normalized weights of i-th particle, respectively.

Step 3, resample particles and calculate the fusion-positioning results. Based on each particle’s normalized weight, the particles with small weights are abandoned and the particles with large weights are copied to replace them. After resampling, the average coordinates of the resampled particles are computed by using Equation (27) and setting as the fusion positionR(k):

R(k)=1n(∑i=1n_xi,∑i=1n_yi)

Step 4, correct the PDR positions. After the three steps above, we should correct the PDR coordinates by using the fusion positionR(k).

Step 5, particles mutation. In this step, we should convert the coordinates of the resampled particles into the binary code by using Algorithm 2. Gene mutation is used to change the coordinate values of the particles and the mutated particles will be passed to the next generation.

To summarize, when the system is running continuously, the above steps will be executed cyclically. At each moment, the system will geneticize the best quality particles to the next moment and the geomagnetic positioning results at the next moment will reconstruct the particle set. Particle mutation and the reconstruction of particle set methods ensure the reliability and the diversity of particles. Figure 9 shows the process of the fusion-positioning system.

3. Experiments and Results

In this section, several experiments are designed to verify the methods mentioned in the above sections. As shown in Figure 10, they are the experimental area scene. The test path connects area A and B and the total length is approximately 130 m. There is a square area which is about 54 square meters between the two areas and the total experimental area is about 360 square meters. We use a Xiaomi 6 mobile phone that has installed our developed data acquisition software to collect data with a sampling rate of 50 HZ. In order to test the performance of the MCF algorithm, the actual azimuth angles from A₁ to A₂, A₂ to B₁, and B₁ to B₂ are measured beforehand, which are 88.1°, 357.17°, and 268.1°, respectively. The size of the particle set is 100 in the genetic-particle filter and all the data processing and simulation work are done on the computers.

3.1. Geomagnetic Fingerprint Database Construction

It is essential for geomagnetic positioning to build a high precision database. In our construction work, we set the reference point every 1.2 m, and then collected the geomagnetic data for 15 s on these points. Five geomagnetic features are extracted by using the method in Section 2.2.1. As for the uncollected area, the linear interpolation method is used to obtain the geomagnetic data. Five geomagnetic maps are generated as illustrated in Figure 11. From these figures, we can find that geomagnetic characteristics of different regions are very obvious, especially Figure 11a which shows the x-axis component is not close to zero in some experimental areas and this is contrary to theory. We analyzed that ferrous materials lead to the anomalous magnetic field. But they can better describe the magnetic field distribution and improve the specificity of the geomagnetic fingerprint. Therefore, x-axis components are extracted and combine the other four features to construct the fingerprint database.

3.2. Geomagnetic Positioning Experiment

In this section, we apply one feature (total intensity of magnetic field), two features (horizontal and total components), three features (three-axis components), four features (three-axis components and horizontal intensity) and five features to geomagnetic positioning. Positioning errors are presented in Figure 12. We also counted the positioning error of 50% and 80% of the test points, and compute the mean positioning error and root mean square error (RMSE) in Table 2.

As shown in Table 2 and Figure 12, one-feature or two-features data can differentiate positions with lower accuracy. When five features data are used, the positioning accuracy and the RMSE are 2.98 m and 4.01 m, which are better than those of using four or fewer features. The positioning accuracy gradually increases with the number of features of data increasing. Based on the experimental results in Section 3.1, five-features data have different distributions in different areas and can better describe the magnetic field at different points. Therefore, we can conclude that five geomagnetic features could improve the geomagnetic positioning accuracy and geomagnetic fingerprint specificity.

3.3. Heading Angle Experiment

In order to evaluate the performance of the Mahony complementary filter (MCF) algorithm, we utilized the mobile phone to collect the heading angles by using the MCF algorithm, electronic compass (EC) and the uncorrected gyroscope on the same path. The estimation errors of these three methods were compared in Table 3 and Figure 13. The mean estimation accuracy and the root mean square error (RMSE) of the MCF algorithm are 7.21 ° and 6.58°, respectively. Compared with the estimation errors of EC and the uncorrected gyroscope, we can conclude that the MCF algorithm can provide more accurate and stable heading angles.

3.4. Fusion-Positioning Experiment In this section, we analyzed the mean positioning accuracy and root mean square error (RMSE) of the pedestrian dead-reckoning (PDR) method, geomagnetic positioning, fusion positioning based on the classic particle filter (PF) which has no gene mutation and particle set reconstruction, and the proposed fusion method. The computational cost (CC) obtained by calculating the average positioning time of different methods is also counted for a comparison.

As shown in Table 4, compared with the PDR method, geomagnetic positioning and fusion positioning by using PF, the mean positioning error of the proposed method is 1.72 m and the RMSE is 1.89 m, respectively. Moreover, 80% of the test points have the positioning error within 2.45 m. But under the same confidence level, the positioning accuracy of geomagnetic positioning, the PDR method and PF algorithm are only 4.51 m, 4.80 m and 3.22 m, respectively. In terms of the computational cost, we can find that our proposed algorithm can give positions within 0.3 s. This is slightly larger than PF, but it is still low because the size of the particle set is just 100. Table 4 shows the positioning accuracy and stability of our proposed algorithm are better improved and the computational cost is low. Moreover, the points positioning errors and the cumulative distribution of positioning errors are presented in Figure 14.

From Figure 14a we can find that the positioning error of all the test points using the proposed fusion-positioning method are within 4 m and more than 70% of test points have the positioning error better than 2.5 m. These are much better than the other three positioning methods. In terms of the points positioning errors, as shown in Figure 14b,c, the positioning error of the PDR method gradually accumulates with time and the error of final point accumulates to more than 6 m. Geomagnetic positioning and the fusion positioning method based on PF are both unstable and several large positioning errors appear in the error curves. Conversely, our proposed fusion-positioning method has neither a large positioning error nor the error accumulation phenomenon. The positioning accuracy and stability are improved.

Based on the above error analysis, we have concluded that our proposed algorithm has better positioning performance. Figure 15 shows the positioning trajectories of different methods, the PDR method can provide more reliable positioning results at first, but the trajectory begins to shift after the first turn. What is worse, the positioning trajectory directly passes through the room after the second turn, which seriously deviates from the actual path. The fusion-positioning method based on PF also shows the instability after the first turn, which suggests that particle degradation affects the positioning stability and accuracy as the time grows. Contrary to the instability and the navigation path drift, the navigation trajectory indicated by our proposed method has no large drift error and is closer to the actual experimental path. There is no doubt that the mutation and reconstruction of the particle set methods improve the performance of the particle filter and ensure the stability and high precision of the fusion-positioning algorithm.

4. Discussion and Conclusions

In this paper, we propose a fusion positioning method based on the genetic-particle filter. In order to optimize the particle degradation problem, we design an optimization mechanism that utilizes the gene mutation method and the method of reconstructing particle set to ensure the reliability and diversity of particles. The experiment results show that the average position estimation error is 1.72 m and the RMSE is 1.89 m, the positioning accuracy and stability are improved. Contrary to the high computational burden caused by increasing the number of particles in papers [26,28,29], the size of the particle set in this paper is 100 and the computational load is lower compared with them. Furthermore, the single-point-based magnetic positioning (SPMP) method is more flexible and easier to be implemented and its computational load is also not very high. Therefore, it is possible for us to implement this fusion positioning algorithm on smartphones. Moreover, five geomagnetic features are extracted in our work and they can improve the fingerprint specificity of a single point and geomagnetic positioning accuracy. Besides, better heading angle estimation performance of the Mahony complementary filter (MCF) algorithm also provides another effective approach for us to calculate the heading angle by using smartphones. All of these results can support that our proposed fusion-positioning method has application value in indoor positioning and navigation.

However, there are also some limitations in our fusion-positioning method. The linear interpolation method can reduce data collection work. Whether there is another effective approach that can reduce more geomagnetic database construction work is still needed to explore in depth. We did not consider the influence of different walking speed on fusion positioning. The effect of different motion states on heading angle estimation is not studied. All of these limitations are our future research work.

Figure 1. The framework of pedestrian dead-reckoning (PDR) technology in this research.

Figure 2.(a) Acceleration curve of the walking process; (b) acceleration curve of one step.

Figure 3. The execution process of the FSM algorithm.

Figure 4.(a) The carrier coordinate system (CCS); (b) roll; (c) pitch; (d) yaw.

Figure 5. Comparison of the untransformed and transformed three-axis geomagnetic data. (a) Untransformed three-axis geomagnetic data; (b) Transformed three-axis geomagnetic data

Figure 6. Gene mutation.

Figure 7. Particle mutation.

Figure 8. Process of genetic-particle filter.

Figure 9. Fusion-positioning algorithm based on genetic-particle filter.

Figure 10.(a) Experimental area; (b) experimental path.

Figure 11. Geomagnetic characteristics interpolation map. (a) X-axis component interpolation map; (b) Y-axis component interpolation map; (c) Z-axis component interpolation map; (d) horizontal component interpolation map; (e) magnetic field strength interpolation map.

Figure 12. Geomagnetic positioning comparison by using different features data.

Figure 13. Heading angles comparison by using different methods (small angles added 360°). (a) Uncorrected gyroscope VS MCF; (b) Electronic compass VS MCF

Figure 14. Positioning cumulative error curve and point error curves of different methods. (a) Positioning cumulative error curves; (b) points error curves of PDR, geomagnetic positioning and the proposed method; (c) points error curves of fusion positioning based on PF and the proposed method.

Figure 15. Positioning trajectories of different methods.

Notation	Description of Notation
acc_diff	Acceleration difference threshold
acc_thr	Entering the FSM algorithm threshold
up_thr	Maximum value judged as a rising state
down_thr	Minimum value judged as a falling state
peak_thr	Peak threshold
D_thr	Max interference value
Pup	Times of acceleration greater thanacc_diff
Pdown	Times of acceleration smaller thanacc_diff
ai	Acceleration value

Features	Min	Max	Mean	RMSE	50% Error	80% Error
One feature	1.16	15.91	8.01	8.74	7.98	10.69
Two features	0.27	16.29	6.62	7.41	7.06	9.54
Three features	0.19	16.27	3.80	4.77	3.15	6.08
Four features	0.02	18.16	3.38	4.66	2.68	4.65
Five features	0.04	12.07	2.98	4.01	2.56	4.51

Algorithm	Min	Max	Mean	RMSE
MCF	0.04	50.78	7.21	6.58
EC	0.05	59.03	11.46	9.81
Gyroscope	0.18	88.05	16.61	13.39

Method	Min/m	Max/m	Mean/m	RMSE/m	50% Error/m	80% Error/m	CC/s
PDR	0.12	6.29	3.14	3.43	2.54	4.80	0.004
Geomagnetic	0.04	12.07	2.98	4.01	2.56	4.51	0.252
PF	0.10	4.83	1.96	2.28	1.84	3.22	0.285
Proposed	0.13	3.42	1.72	1.89	1.75	2.45	0.296

Author Contributions

Conceptualization, M.S. (Meng Sun); Formal analysis, M.S. (Meng Sun) and S.X.; Funding acquisition, Y.W.; Methodology, H.C.; Project administration, Y.W.; Data curation, H.C., M.S. (Minghao Si); Writing-original draft preparation, M.S. (Meng Sun). All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China [grant number 2016YFB0502102], in part by the State Key Laboratory of Satellite Navigation System and Equipment Technology [grant number CEPNT-2018KF-03].

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the writing of the manuscript; and in the decision to publish the results.

References

1. Davidson, P.; Piché, R. A survey of selected indoor positioning methods for smartphones. IEEE Commun. Surv. Tutor. 2017, 19, 1347-1370.

2. Li, X.; Zhang, X.; Ren, X.; Fritsche, M.; Wickert, J.; Schuh, H. Precise positioning with current multi-constellation global navigation satellite systems: GPS, GLONASS, Galileo and BeiDou. Sci. Rep. 2015, 5, 8328.

3. Lin, K.; Chen, M.; Deng, J.; Hassan, M.M.; Fortino, G. Enhanced fingerprinting and trajectory prediction for IoT localization in smart buildings. IEEE Trans. Autom. Sci. Eng. 2016, 13, 1294-1307.

4. He, S.; Chan, S.-H.G. Wi-Fi fingerprint-based indoor positioning: Recent advances and comparisons. IEEE Commun. Surv. Tutor. 2016, 18, 466-490.

5. Yoon, P.K.; Zihajehzadeh, S.; Kang, B.-S.; Park, E.J. Adaptive Kalman filter for indoor localization using Bluetooth Low Energy and inertial measurement unit. In Proceedings of the 2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), Milan, Italy, 25-29 August 2015; pp. 825-828.

6. De Blasio, G.; Quesada-Arencibia, A.; García, C.R.; Rodríguez-Rodríguez, J.C.; Moreno-Díaz, R. A protocol-channel-based indoor positioning performance study for Bluetooth low energy. IEEE Access 2018, 6, 33440-33450.

7. Schroeer, G. A Real-Time UWB Multi-Channel Indoor Positioning System for Industrial Scenarios. In Proceedings of the 2018 International Conference on Indoor Positioning and Indoor Navigation (IPIN), Danvers, MA, USA, 24-29 September 2018; pp. 1-5.

8. Chen, P.; Ye, K.; Chen, X. A UWB/Improved PDR Integration Algorithm Applied to Dynamic Indoor Positioning for Pedestrians. Sensors 2017, 17, 2065.

9. Bekkali, A.; Sanson, H.; Matsumoto, M. RFID indoor positioning based on probabilistic RFID map and Kalman filtering. In Proceedings of the Third IEEE International Conference on Wireless and Mobile Computing, Networking and Communications (WiMob 2007), New York, NY, USA, 8-10 October 2007; p. 21.

10. Yucel, H.; Edizkan, R.; Ozkir, T.; Yazici, A. Development of indoor positioning system with ultrasonic and infrared signals. In Proceedings of the 2012 International Symposium on Innovations in Intelligent Systems and Applications, Trabzon, Turkey, 2-4 July 2012; pp. 1-4.

11. Huang, C.-N.; Chan, C.-T. ZigBee-based indoor location system by k-nearest neighbor algorithm with weighted RSSI. Procedia Comput. Sci. 2011, 5, 58-65.

12. Davidson, P.; Takala, J. Algorithm for pedestrian navigation combining IMU measurements and gait models. Gyroscopy Navig. 2013, 4, 79-84.

13. Zheng, L.; Zhou, W.; Tang, W.; Zheng, X.; Peng, A.; Zheng, H. A 3D indoor positioning system based on low-cost MEMS sensors. Simul. Model. Pract. Theory 2016, 65, 45-56.

14. Siiksakulchai, S.; Thongchai, S.; Wilkes, D.M.; Kawamura, K. Mobile robot localization using an electronic compass for corridor environment. In Proceedings of the 2000 IEEE International Conference on Systems, Man, and Cybernetics, Nashville, TN, USA, 8-11 October 2000; pp. 3354-3359.

15. Li, B.; Gallagher, T.; Dempster, A.G.; Rizos, C. How feasible is the use of magnetic field alone for indoor positioning? In Proceedings of the 2012 International Conference on Indoor Positioning and Indoor Navigation (IPIN), Nashville, TN, USA, 8-10 October 2000; pp. 1-9.

16. He, S.; Shin, K.G. Geomagnetism for smartphone-based indoor localization: Challenges, advances, and comparisons. ACM Comput. Surv. (CSUR) 2018, 50, 97.

17. Senin, P. Dynamic time warping algorithm review. Inf. Comput. Sci. Dep. Univ. Hawaii Manoa Honol. USA 2008, 855, 40.

18. Subbu, K.P.; Gozick, B.; Dantu, R. LocateMe: Magnetic-fields-based indoor localization using smartphones. ACM Trans. Intell. Syst. Technol. (TIST) 2013, 4, 73.

19. Subbu, K.P.; Gozick, B.; Dantu, R. Indoor localization through dynamic time warping. In Proceedings of the 2011 IEEE International Conference on Systems, Man, and Cybernetics, Anchorage, AK, USA, 9-12 October 2011; pp. 1639-1644.

20. Ashraf, I.; Hur, S.; Park, Y. mPILOT-Magnetic Field Strength Based Pedestrian Indoor Localization. Sensors 2018, 18, 2283.

21. Ashraf, I.; Hur, S.; Shafiq, M.; Park, Y. Floor Identification Using Magnetic Field Data with Smartphone Sensors. Sensors 2019, 19, 2538.

22. Hui, L.; Luo, H.; Fang, Z.; Li, X. TACO: A Traceback Algorithm Based on Ant Colony Optimization for Geomagnetic Positioning; Springer: Berlin/Heidelberg, Germany, 2014; pp. 208-224.

23. Shu, Y.; Shin, K.G.; He, T.; Chen, J. Last-mile navigation using smartphones. In Proceedings of the 21st Annual International Conference on Mobile Computing and Networking, Paris, France, 7-15 September 2015; pp. 512-524.

24. Le Grand, E.; Thrun, S. 3-axis magnetic field mapping and fusion for indoor localization. In Proceedings of the 2012 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI), Hamburg, Germany, 13-15 September 2012; pp. 358-364.

25. Kim, B.; Kong, S.H. A Novel Indoor Positioning Technique Using Magnetic Fingerprint Difference. IEEE Trans. Instrum. Meas. 2016, 65, 2035-2045.

26. Huang, H.; Li, W.; Luo, D.A.; Qiu, D.W.; Gao, Y. An Improved Particle Filter Algorithm for Geomagnetic Indoor Positioning. J. Sens. 2018, 2018, 5989678.

27. Benesty, J.; Chen, J.; Huang, Y.; Cohen, I. Pearson correlation coefficient. In Noise Reduction in Speech Processing; Springer: Berlin, Germany, 2009; pp. 1-4.

28. Zheng, M.; Zhao, Y.; He, Y.; Xu, C.-Z. Sensitivity-based adaptive particle filter for geomagnetic indoor localization. In Proceedings of the 2017 IEEE/CIC International Conference on Communications in China (ICCC), Qingdao, China, 22-27 October 2017; pp. 1-6.

29. Shu, Y.; Bo, C.; Shen, G.; Zhao, C.; Li, L.; Zhao, F. Magicol: Indoor localization using pervasive magnetic field and opportunistic WiFi sensing. IEEE J. Sel. Areas Commun. 2015, 33, 1443-1457.

30. Srinivas, M.; Patnaik, L.M. Genetic algorithms: A survey. Computer 2002, 27, 17-26.

31. Mahony, R.; Hamel, T.; Pflimlin, J.-M. Nonlinear Complementary Filters on the Special Orthogonal Group. IEEE Trans. Autom. Control 2008, 53, 1203-1218.

32. Butcher, J.C.; Butcher, J. The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods; Wiley-Interscience: New York, NY, USA, 1987; Volume 512.

33. Ahn, D.; Park, J.-S.; Kim, C.-S.; Kim, J.; Qian, Y.; Itoh, T. A design of the low-pass filter using the novel microstrip defected ground structure. IEEE Trans. Microw. Theory Tech. 2001, 49, 86-93.

34. Keselbrener, L.; Keselbrener, M.; Akselrod, S. Nonlinear high pass filter for R-wave detection in ECG signal. Med. Eng. Phys. 1997, 19, 481-484.

35. Zhang, P.; Li, F. A new adaptive weighted mean filter for removing salt-and-pepper noise. IEEE Signal Process. Lett. 2014, 21, 1280-1283.

36. Brajdic, A.; Harle, R. Walk detection and step counting on unconstrained smartphones. In Proceedings of the 2013 ACM International Joint Conference on Pervasive and Ubiquitous Computing, Zurich, Switzerland, 8-12 September 2013; pp. 225-234.

37. Zhuang, Y.; Lan, H.; Li, Y.; El-Sheimy, N. PDR/INS/WiFi integration based on handheld devices for indoor pedestrian navigation. Micromachines 2015, 6, 793-812.

38. Ladetto, Q. On foot navigation: Continuous step calibration using both complementary recursive prediction and adaptive Kalman filtering. In Proceedings of the ION GPS, Salt Lake City, UT, USA, 19-22 September 2000; pp. 1735-1740.

39. Vildjiounaite, E.; Malm, E.-J.; Kaartinen, J.; Alahuhta, P. Location estimation indoors by means of small computing power devices, accelerometers, magnetic sensors, and map knowledge. In Pervasive Computing; Springer: Berlin, Germany, 2002; pp. 211-224.

40. Tian, Q.; Salcic, Z.; Kevin, I.; Wang, K.; Pan, Y. A multi-mode dead reckoning system for pedestrian tracking using smartphones. IEEE Sens. J. 2015, 16, 2079-2093.

41. Weinberg, H. Using the ADXL202 in pedometer and personal navigation applications. Analog Devices AN-602 Appl. Note 2002, 2, 1-6.

42. Blu, T.; Thévenaz, P.; Unser, M. Linear Interpolation Revitalized. IEEE Trans. Image Process. A Publ. IEEE Signal Process. Soc. 2004, 13, 710-719.

43. Deng, Z.; Zhu, X.; Cheng, D.; Zong, M.; Zhang, S. Efficient kNN classification algorithm for big data. Neurocomputing 2016, 195, 143-148.

44. Ma, Y.; Dou, Z.; Jiang, Q.; Hou, Z. Basmag: An optimized HMM-based localization system using backward sequences matching algorithm exploiting geomagnetic information. IEEE Sens. J. 2016, 16, 7472-7482.

45. Arulampalam, M.S.; Maskell, S.; Gordon, N.; Clapp, T. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal Process. 2002, 50, 174-188.

46. Gay, D.M. Correctly rounded binary-decimal and decimal-binary conversions. Numer. Anal. Manuscr. 90-10. 1990. Available online: https://ampl.com/REFS/rounding.pdf (accessed on 29 November 2019).