Integrity Monitoring for Bluetooth Low Energy Beacons RSSI Based Indoor Positioning

Indoor wireless positioning using Bluetooth Low Energy (BLE) beacons have attracted considerable attention from industry and academia given the many advantages of this technology such as low power consumption, low cost, easy layout, high availability, and high precision. However, the indoor positioning accuracy always suffers from non-line of sight (NLOS) propagation, stemming from the frequently occurring instances of reflection, refraction, or scattering of BLE radio signals due to the complexity of indoor environments. This article proposes an integrity monitoring (IM) algorithm to detect and eliminate two gross errors simultaneously to solve the adverse effects caused by the NLOS. The logarithmic attenuation model based on the received signal strength indication (RSSI) of BLE realizes positioning by combining trilateration and Least Squares Based on the Taylor expansion (LSBT). Furthermore, the IM based on hypothesis testing is employed to improve the positioning quality andthe users will be alerted in time to avoid risk from positioning accuracy no longer meet user’s requirement. The performance of the proposed IM algorithm has been extensively tested by conducting simulation and field experiments. The experimental results show that the IM algorithm significantly improved BLE positioning accuracy as well as the robustness of the positioning system. The 90% average error (1.9143m) in seven groups of single point experiments was reduced by 34.48% over the 90% average error (2.9143m) of the LSBT method after performing IM, and the maximum error during continuous positioning did not exceed 3m after performing IM, which were better than only using LSBT positioning.


I. INTRODUCTION
With the rapid development in the performance of the Global Navigation Satellite System (GNSS) and indoor positioning technology, location-based services (LBS) have progressed fast in last decades and have become an indispensable part of everyday life. Positioning and navigation systems can give users an accurate description of their current position, timing, and posture in a space-time coordinate system, and thus meet user demands in diverse application cases such as transportation, construction engineering, smart city plan-The associate editor coordinating the review of this manuscript and approving it for publication was Kegen Yu . ning, and emergency response. However, a positioning and navigation system is not always in a stable state. In practice, the positioning sensors have faults or the signals are disturbed during the propagation process, they may cause a serious deviation in the positioning results. The faults defined as the gross errors in the observations during positioning process, while the observations without faults only contain random noise. If the positioning and navigation system fails to exclude faults or to alert users in time when faults appear in a positioning system or when the positioning accuracy no longer matches the application requirement, it will provide hazardous misleading information to those users. Therefore, the integrity of the positioning and navigation system has VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ become an indispensable indicator when evaluating system performance.
Integrity is one of the four navigation performance indicators, i.e., accuracy, integrity, continuity, and availability [1]. It detects and identifies faults or abnormality that affect positioning accuracy and to inform the user that the service is unavailable. Integrity is a confidence measure referring the system ability to provide timely notifications or to terminate a signal when the positioning and navigation system is unavailable, due to faults or when positioning errors exceed the alert limit [2].
Various indoor applications, e.g., firefighting, peacekeeping missions, emergency safety, disaster relief, and mobile health services and intelligent robot applications, require the positioning system with high integrity. Unfortunately, indoor environments are more complex and varied than in outdoor environments. Radio waves are reflected, refracted, or scattered frequently due to the blocked by obstacles, which change the signal propagation path to the receiver forming a non-line of sight (NLOS) propagation and multipath effects. The NLOS propagation and multipath effects will result in a large deviation in the positioning results that can seriously affect positioning accuracy [3]. Moreover, the indoor spatial layout and topology are susceptible to human influences such as changes in sound, light, and the electromagnetic environment in indoor spaces, for positioning methods based on feature matching, the positioning results will be greatly affected [4].
It is challenging to realize high precision, highly robust indoor positioning system. Nevertheless, integrity monitoring is a desirable option to improve positioning system performance. But most of the existing researches about integrity monitoring algorithm of indoor positioning have some limitations [5]- [9]: Only single fault (gross error) can be identified and eliminated during the positioning process (that is, these integrity monitoring algorithms cannot identify and eliminate multiple gross errors during the positioning process); furthermore, these integrity monitoring algorithms all lack an alert mechanism for cases in which positioning accuracy completely dissatisfies user requirements. However, since indoor environments are complicated by reflection, shadowing and multipath, multiple gross errors occur frequently during indoor positioning. Besides, if an application system only outputs information but lacks a reliable description of that information, the service is incomplete essentially. Therefore, when the positioning system is unavailable or the positioning accuracy does not meet user requirements, the positioning system should alert the user in time.
Currently, there are many positioning sources used by indoor positioning technology: WiFi [10], Pedestrian Dead Reckoning (PDR) [11], Bluetooth Low Energy (BLE) [12], ZigBee [13], RFID [14], Ultra Wide Band (UWB) [15], Ultrasound [16], Infrared Signal [17], the Computer Vision [18], the Magnetic field [19], Optical localization [20] and so on. Among those indoor positioning technologies, BLE positioning has many advantages such as low power consumption, low cost, easy layout, high availability, and high precision. The low-energy and low-cost Bluetooth devices and the performance of mobile intelligent terminals have experienced a fast progress in the last decade. In 2010, the release of BLE protocol (Bluetooth 4.0 standard) effectively reduced energy consumption; In 2016, Bluetooth 5.0 standard was released, adding indoor positioning assistance function, positioning distance up to 200m and with less power consumption [21].
Bluetooth has higher precision and saves more energy than WiFi. In contrast to ZigBee and RFID, Bluetooth is easier to deploy. As compared to UWB, Ultrasound, Infrared signals, Computer vision, Magnetic field, and Optical localization, the Bluetooth costs less and more easily implemented [22]. Moreover, all smartphones and wearable devices have a Bluetooth module. Because of these advantages, Bluetooth has great potential for indoor positioning [23]- [30].
At present, indoor positioning algorithms based on Bluetooth can be mainly divided into two categories [21]: geometric methods and methods based on Received Signal Strength Indication (RSSI). This article is based on the RSSI methods. The geometric methods like Time of Arrival (TOA) [31], [32], Angle of Arrival (AOA) [33], and Time Difference of Arrival (TDOA) requiring specialized hardware, which is costlier than RSSI methods [34]. The RSSI based positioning methods commonly used can be divided into range-free methods and range-based methods [35], a classic range-free method based on RSSI is the fingerprinting positioning method as in [26], [28]. This method has high accuracy, but needs to collect many fingerprint information offline, which will require a long training period. It can be realized in a small indoor environment and is not easy to realize in a large indoor environment. Besides, the fingerprint information needs to update periodically to guarantee positioning accuracy. In contrast, the range method based on RSSI does not need to collect fingerprint information offline and uses the wireless signal propagation model to convert the Bluetooth signal strength into the distance. The range method based on RSSI is low-cost and easy to implement. In this article, we propose an algorithm for smartphone-based indoor positioning using RSS signals from BLE beacons. The algorithm estimates positions using the distances between the target and BLE beacons, are made by combining trilateration and the Least Squares Based on the Taylor expansion (LSBT).
To further improve accuracy and robustness, we propose a new integrity monitoring (IM) algorithm for RSSI-based indoor positioning using BLE Beacons. In our research, the major innovation is the IM alert mechanism that alerts users timely for reduce the risk from positioning accuracy no longer meet user's requirements. We collected the RSSI fitting the Bluetooth transmission and the distance parameter of the logarithmic attenuation model, arriving at a suitable Bluetooth ranging model that combines trilateration and LSBT to realize indoor positioning. This integrity monitoring algorithm based on hypothesis testing to simultaneously identify and remove two gross errors during the positioning process. The integrity monitoring constituted by the parity vector method and the Maximum Likelihood Estimation Gross Error Test (MLEGET) method. We conducted simulations and a series of field experiments including single point positioning and continuous positioning to test and verify the proposed algorithm. The experimental results show that the proposed IM algorithm can effectively improve the positioning accuracy and system robustness. The users will be alarmed in time to avoid risk from the positioning accuracy does not meet the user requirements due to positioning fault.
The rest part of this article is arranged as follows. Section II reviews the related work about integrity monitoring. Section III states the methodology for BLE RSSI ranging model, indoor positioning algorithm based on ranging and integrity monitoring algorithm. Section IV and Section V construct the simulation analysis and experimental evaluation which verifies the correctness of positioning theory. Section VI is discussion and conclusions.

II. RELATED WORK
In an outdoor environment, different GNSS integrity indicators have been developed to meet various application requirements. The International Civil Aviation Organization(ICAO) sets a multi-level integrity limit for aircraft from remote ocean areas to precise approach phases. The vertical alarm limit is 40 meters during a precision approach and the horizontal alarm limit is 15 meters, with an alarm time of fewer than six seconds [36]; The Galileo system sets requirements for life safety and public control services based on the integrity risk value [37]. The Galileo system judges the integrity content by monitoring the Signal in Space Error (SISE), Signal in Space Accuracy (SISA), Signal in Space Monitoring Accuracy (SISMA), and sends the monitoring results to users. Besides, the Galileo system provides five types of services face to global: Open Services, Safety of Life, Commercial Service, Public Regulation Services, Search and Rescue. These services are all related to the integrity monitoring of the Galileo system.
Among research on integrity monitoring related to GNSS applications, the Receiver Autonomous Integrity Monitoring (RAIM) algorithm is widely used. This integrity monitoring algorithm on the client-side has many advantages, such as independence from external devices, low cost, and easy to implement. The classical integrity monitoring algorithms include pseudo-range comparison, parity vector and leastsquares residuals [38].
The indoor environment is usually much complex than the outdoors with severe multipath propagation, and wireless signal interruption, which makes the GNSS integrity monitoring algorithm not fully applicable in indoors scenarios. However, the concept of integrity as proposed in the context of GNSS is still applicable to the performance evaluation of indoor positioning systems.
Unlike GNSS, research in indoor positioning currently main focus on positioning accuracy rather than advanced system performance such as integrity monitoring. Only few researches on indoor positioning integrity monitoring algorithms have been reported. Yerubandi et al. [5] proposed an integrity monitoring algorithm to assess the continuous availability of WLAN access points (APs), which can identify and isolate rogue APs to improve positioning accuracy, by taking account of integrity monitoring, the performance and robustness of the indoor positioning system were improved. Similarly, in [6], a mobile device WLAN indoor positioning algorithm based on integrity monitoring was proposed, this integrity monitoring algorithm utilizes redundancy at the APs and isolates APs with the damaged features, which improves the robustness of the system. Ascher et al. [7] investigated a UWB/INS positioning system and proposed an integrity monitoring algorithm based on innovation of extended Kalman filter, which can effectively detect and isolate the TDOA observation outliers to further improve positioning accuracy. Wang et al. [8] proposed a method of using the short-time reliability of PDR to aid in visual integrity monitoring and to reduce positioning error. Zhao et al. [9] aimed at indoor high-precision engineering application proposed carrier phase-based integrity monitoring (CRAIM) algorithm to identify and exclude potential faults of the pseudolites, the CRAIM ensured the accuracy and reliability of positioning results and achieved accuracies at the centimeter level for dynamic experiments and millimeter levels for static ones. In [39], specific thresholds of integrity monitoring were proposed to the prevention of theft and criminal surveillance, etc., including the alarm time (TTA, Time to Alarm) and the Horizontal Alarm Limit (HAL). A partial list of specific indicators is shown in Table 1. However, the study did not carry out specific integrity monitoring algorithms.
In [5]- [9], the proposed integrity monitoring algorithms only focus on a single fault (single gross error), and all lack an alert mechanism for cases when positioning accuracy completely dissatisfies user requirements.

III. METHODOLOGY
This section first presents the BLE RSSI ranging model for distance estimation, which is followed by trilateration and LSBT for positioning estimation. Next, the details of the IM algorithm are described to remove gross errors and improve the performance of an indoor positioning system. Finally, the availability of the IM algorithm is evaluated to determine whether to alert users in time to avoid risk from positioning accuracy no longer meet user's requirements.

A. BLE RSSI RANGING MODEL
Numerous theoretical derivation and empirical formulas [12], [35], [40] show that there is a definite relationship between the RSSI and range. In this article, we analyze the relationship between BLE RSSI and transmission distance based on the logarithmic attenuation model in [12], [35], [40] as shown in (1).
where P (d) represents the RSS at the distance d between the transmitting end and the receiving end; P d 0 represents the RSS at the reference distance d 0 ; n is the path-loss exponent; ε represents a Gaussian random variable, with zero means, caused by shadow fading [40].

B. INDOOR POSITIONING ALGORITHM BASED ON RANGING
In our research, the initial value of the BLE positioning point to be measured is calculated by the trilateration method, and the BLE positioning is realized by the LSBT. The position of a point can be calculated using the trilateration method [41]. Given the known coordinates of at least three non-linear BLE beacons and the corresponding distances between the BLE beacons and the point, as shown in Fig.1. The coordinates (x, y) of a point to be measured can be solved by the following expression: where (x 1 , y 1 ) , (x 2 , y 2 ) . . . (x n , y n ) are the coordinates of BLE beacons, d 1 , d 2 . . . d n denote the distances between the beacons and the target node. linearizing the Equation (2) based on the Taylor formula (i.e. using LSBT) at the initial value x 0 , y 0 , one can get: where ( x, y) is the coordinate estimation, d is the distance observation vector, d 0 is the distance vector calculated based on the initial coordinate value of the point to be measured, B denotes the design matrix, ε represents the random error vector, where ε ∼ N (0, σ 2 I n×n ), X represents the estimated coordinate value of the point to be measure, X 0 represents the initial coordinate value, the initial coordinate value x 0 , y 0 of the point to be measured can be obtained by the trilateration method given three non-linear BLE beacons and corresponding distances.
The estimated value x of correction x can be calculated based on the least-squares residual squared sum minimum criterion as: bring the x into equation (3) to calculate residual vector v as: [42]. The estimated coordinate value X can be obtained by adding the x and the initial value X 0 , and then take the obtained coordinate value as the new initial value for iterative calculation. When the absolute value of x is less than a certain threshold, iteration is stopped and the final solution is obtained.

C. INTEGRITY MONITORING ALGORITHM
Since indoor environments are very complicated, which will result in BLE radio signals form NLOS propagation and multipath effects. Consequently, multiple gross errors occur frequently during indoor positioning lead to the positioning accuracy relatively low realized by the trilateration and the LSBT. To further improve positioning accuracy and robustness of the positioning system, a new integrity monitoring algorithm based on the parity vector for single fault (gross error) detection is proposed in this section. This algorithm can simultaneously identify and eliminate two gross errors. The algorithm can also alert the user in time when the positioning accuracy does not meet the user requirements.

1) PARITY VECTOR
The parity vector method is a chi-square test based on a single gross error hypothesis, as proposed by Sturza in 1988 [42]. The process is given below.
Assume that all observations have the same observation accuracy, and the observation equation is as shown in Equation (3). QR decomposition of the design matrix B is obtained: where n represents the aggregate of BLE beacons observations participating in the positioning at the current time, and u is the necessary number of BLE beacons observations for calculating the coordinate estimation value, Q is an identity orthogonal matrix of order n × n, and R is (n × u) an upper triangular matrix, and Q p represents the parity transformation matrix, whose rows are mutually orthogonal, unitary in magnitude, and orthogonal to the columns of B, and Q p B = 0 and Q T p Q p = S. The two equations are derived in Appendix. Projecting the observation error into the parity space, we can get the parity vector p as: since parity vector p reflects the observation error information, test statistics can be constructed based on it to carry out fault detection and identification. The fault test statistics using the parity vector construction is as follows: where σ is the standard deviation of the observations. During the positioning process, when the observation error ε only contains an random error, namely ε ∼ N (0, σ 2 I n×n ), the test statistics T obeys the center χ 2 distribution with a degree of freedom (DOF) of n − u [42]: When the observation error ε i contains the gross error ∇b i , the test statistics T obeys the non-central χ 2 distribution with a DOF of n − u: where λ is a non-central parameter. According to the idea of hypothesis testing, a fault detection threshold T C can be calculated from the cumulative probability density function of the center χ 2 distribution when the probability of false alarm P FA and DOF are given. The non-central parameter λ can be obtained from the cumulative probability density function of the non-central χ 2 distribution when the probability of missed detection P MD and DOF are given, as follows: Taking the chi-square distribution with a DOF of 2 as an example, the relationship between P FA , P MD , and T C is as shown in Fig. 2. According to the Fig. 2, the P FA and P MD are correlated, i.e., the lager the P FA , the smaller the P MD and the T C are. According to the relationship between the non-center parameter λ and the random variable expectation [43], the non-central parameter λ also can be expressed as follows: where µ is the expectation of the random variable (i.e. the test statistics T ), S matrix is an idempotent matrix, as can be seen from Section III.C.1). The expected E (v) of the residual v is zero when the observation error vector ε contains only random errors; The expectation of v is E (v) = S∇b when the observed error vector ε contains both the random error and the gross error vector ∇b. Therefore, when the i-th row element ∇b i in the gross error vector ∇b is not zero, the relationship between the λ and the ∇b can be expressed: where σ denotes the standard deviation of observations, S ii is the i-th diagonal element of S. The lower bound of the gross error detection of the parity vector method under the single alternative hypothesis can be written as the formula (19) through the equation (18): According to the formula (19), the lower bound of the gross error detection is affected by three factors: 1) The accuracy of observations represented by the standard deviation σ of observations. 2) The geometry of sensor networks in the VOLUME 8, 2020 current time positioning environment, which affects positioning performance and 3) The probability of false alarm P FA and the probability of missed detection P MD , which together determine the non-center λ.
The fault test statistic T is calculated during the positioning process, the gross error is considered to exist when the T is greater than the fault detected threshold T C . Further, it is necessary to identify which observations of the gross error occurred. Considering that the number of observations n is 5, the necessary number of observations u is 3, the impact of gross error on the parity vector p can be obtained according to the equation (11): when the gross error occurs in the observation d i , the projection length of the parity vector p on the i-th column of the parity matrix Q p should be the maximum, so that the observation including the gross error can be identified and then eliminated, and the coordinate values of the point to be measured can be estimated again using the LSBT method.

2) MAXIMUM LIKELIHOOD ESTIMATION GROSS ERROR TEST
The parity vector is a kind of chi-square test method for single gross error detection. If the gross error is regarded as a result of a small probability event in the positioning process, the event within two gross errors simultaneously can satisfy most of the scenarios. Based on this, a detection method for simultaneously detecting two gross errors is proposed in this article.
When there are two gross errors at the same time in the positioning process, the test statistic T in Section III.C.1) also obeys the non-central χ 2 distribution. When there are two gross errors, the relationship between the non-central parameter λ and the gross error ∇b is: where S ii and S jj are the i-th and j-th diagonal element of S, and i = j. When the number of gross errors is unknown, the parity vector method is employed to detect and identify single gross error firstly. if the test statistic T greater than test threshold T C , the observation value d i can be regarded as containing the gross error, the value d i is removed and the LSBT is performed again, and then the T is calculated once again. If the T is still greater than T C , which means the first time detect and identify for gross error is incorrect, that is, there is existing an incorrect gross error identification so that not eliminate the real gross error. It is considered as probably containing two gross errors. Then the parity vector p can be written as: where i = 1, 2, . . . , n; j = 1, 2, . . . , n. The random error ε∼N (0,σ 2 I n×n ). Therefore, according to the idea of maximum likelihood estimation, the likelihood equation is: is the smallest, the likelihood function value is the largest. Therefore, the estimated value ∇b of ∇b can be obtained as: bring ∇b into equation (22) to calculate the estimated value p of the parity vector p as: the likelihood equation (23) can be rewritten as: (27) calculate all the gross error combinations ∇b i , ∇b j , and the number of combinations is C 2 n . when p T p = max p T p , d i and d j are observations with gross errors, and eliminated, The LSBT is performed again. In particular, when the number of gross errors is one, this criterion is equivalent to the parity vector method.
The proposed IM algorithm is only available when certain conditions are met. Section III.D evaluate its availability judgment.

D. INTEGRITY MONITORING AVAILABILITY JUDGEMENT
The IM availability assessment is divided into two steps: the first step determines whether the total number n of BLE beacons meet the requirement, and the second step determines whether the IM algorithm is valid according to the Horizontal Alarm Limit(HAL). The parity vector method can be used for fault detection only when the observations are redundant. In this study, the indoor positioning focus on the two-dimensional planar positioning, so the number of necessary BLE beacons observations is three. Therefore, to make observations redundant, the number of available BLE beacons is at least four. In the first step, as long as the necessary number of BLE beacons participating in positioning is u = 3, and nu = 1, it is possible to detect whether a gross error occurs in the positioning result. when nu > 1, the positioning results of the gross error can be identified. Therefore, when the n > u+1, the second step can be executed; otherwise, the IM is directly considered as unavailable. The specific process is shown in Fig.3. The Fig.3 depicts the detailed IM availability evaluation process. First, the observations redundancy is evaluated. The Horizontal Protection Level (HPL) is calculate and compared to the HAL. If HPL is greater than HAL, the IM unavailability, the positioning system should be alerted to the user in time for reduce the risk from positioning accuracy no longer meet user requirements. Otherwise, the IM is available. Finally, by constructing test statistics to determine whether there is a fault (gross error).
For each currently observed BLE beacons during positioning process, there is a corresponding characteristic slope line, which is the function of the design matrix B. When the user moves within the positioning area, the BLE beacons characteristic slope SLOPE (i) can be written as: For a given horizontal radial error (HRE), the BLE beacon with the maximum characteristic slope SLOPE max , which has the minimum test statistic, and also it is the most difficult to be detected. The HRE can be calculated as following [44]: Fig.4 describes the data cloud using the oval-shaped scatter, whose center that represents the deterministic bias plot without noise will be on the line of SLOPE max if the bias occurs on the BLE beacon with the maximum characteristic slope. The brown-yellow dotted line on the right side in Fig.4 corresponds to the minimum detected bias pbias in the case of the probability of false alarm and probability of missed detection in the parity space. The pbias is related to the Gaussian distribution observation error ε ∼ N (0, σ 2 I n×n ), the geometry of sensor networks, and the probability of false alarm and probability of missed detection associated with the application scenario, the pbias can be written as follows [44]: The HRE corresponding to the intersection point of SLOPE max line and the pbias is the HPL (red dotted line in Fig.4). The calculation formula for HPL can be written as: As shown in Fig.4, the SLOPE (i) of each BLE beacons correspond to its HPL, but the line with SLOPE max is the maximum HPL of all BLE beacons at the current time.
The HAL is given according to the actual application requirements, therefore, whether IM is available or not can be judged according to the following criteria: when the HPL is greater than or equal to the HAL, the IM algorithm is unavailable and the user should be alarmed in time.

IV. SIMULATION ANALYSIS
In this section, The MLEGET was used as the IM algorithm in the simulation experiment. First, the simulation analyzed the impact of Horizontal Dilution of Precision (HDOP) [45] on HPL. Then, the impact of HDOP on the gross error detection performance of the algorithm proposed was analyzed. Finally, we verified the effectiveness of the gross error detection performance of the IM algorithm by setting different gross error magnitude. The simulation environment was a rectangular indoor area with 100m × 50m. The simulated parameters were used from the parameters for the integrity of a precise approach, which set by the ICAO. False alarm rate P FA = 0.333×10 −6 , missed detection rate P MD = 0.001, observation error, ε ∼ N 0,σ 2 I n×n , σ = 0.5m.

A. THE IMPACT OF HDOP ON HPL
According to the analysis of the HPL presented in Section III.D., the positioning performance is affected by the geometry of sensor networks. The plane precision factor can be expressed by the HDOP [45]: where tr refers to the track of the matrix. In theory, the product of HDOP and the ranging accuracy σ express the estimation accuracy of the plane coordinate of the target node [45]. A smaller HDOP value means a more optimal the geometry of sensor networks for positioning. As the number of BLE beacons n increases, the average HDOP in the measurement range decreases, the BLE beacons networks geometry is more beneficial for obtaining higher horizontal positioning accuracy [45]. The HPL values were calculated for the simulation environment with n = 4, 6, 8, 10. The results are shown in Fig.5 and Table 2.  As can be seen from Fig.5, when BLE beacons are evenly evenly deployed, the HPL values gradually increased from the middle area to the surroundings. According to table 2, the maximum value of HPL was approximately twice the average value of HPL. The average value and maximum value of the HPL gradually decreased as the number of visible BLE beacons increased in the positioning area. IM needs to be 100% available to simulate IM gross error detection performance. Therefore, we selected a moderate number of BLE beacons n = 6 for the next simulation experiments.

B. THE IMPACT OF HDOP ON GROSS ERRORS DETECTION PERFORMANCE
The lower bound for gross error detection and identification are related to the BLE beacons network geometry, according to the Section III.C.1). The relationship between the HDOP and gross error detection performance was analyzed through a simulation in this section. Two indicators were selected here to evaluate gross error detection performance: the gross error detection rate and gross error identification rate. The gross error detection rate refers to the probability that the positioning system detects all gross errors when IM is available. The gross error identification rate is the probability of correctly identifying all gross errors during the positioning process when gross errors are detected. The HDOP values distribution when the numbers n of BLE beacons equal 4, 6, 8 and 10 are shown in Fig.6. In this simulation experiment, the number of BLE beacons was set to six, thus, the HDOP distribution values of the experimental area as shown in Fig.6 (b). As can be seen from Fig.6 (b), when the target node moved along the center to the four corners, the HDOP value gradually increased. Therefore, from the center point to the BLE beacons in the lower-left corner, 50 sample points were systematically selected; 50 points were collected on the diagonal from the center of the rectangle to the lower-left corner, and each point was simulated 1000 times. The gross error was set to 5 meters. The probability of gross error appearing in each observation was the same. A single gross error and two gross errors solutions were selected for this simulation. The relationship between HDOP and gross errors detection rate and gross errors identification rate are shown in Fig.7 and Fig.8 respectively.
According to Fig.7 and Fig.8, and the simulation experimental statistical results, it found that the minimum of single gross error detection rate was 97.4%, and the average value of single gross error detection rate was 98.9%. The gross error detection rate decreased as HDOP increased and the  gross error identification rate is close to 100%; The minimum of the two gross errors detection rate was 77.7%, and the average value of the two gross errors detection rate was 85.9%. The two gross error detection rate also decreased as HDOP increased; The gross error identification rate increases with HDOP increased, and remained stable, with an average identification rate of 96.8%. Therefore, as HDOP increased, that is, the observation geometric conditions became deteriorated, the gross error detection rate gradually decreased, and the identification rate of gross error gradually increased and remained stable.

C. THE IMPACT OF THE DIFFERENT MAGNITUDE OF A GROSS ERROR ON GROSS ERRORS DETECTION PERFORMANCE
This section assesses the impact of different magnitudes of gross error on the detection performance, given the target node was fixed. The coordinate of target node was set to (30,20) in Fig.6 (b), the magnitude of gross error varied from 1 to 10 m, and the interval of the magnitude was 0.2 m. The other simulation variables remained the same as those in Section IV.B, each magnitude of the gross error was simulated 10,000 times. The simulation results are shown in Fig.9 and Fig.10.
According to the Fig.9, the average identification rate of IM for a single gross error was 99.76%, and the minimum identification rate was 96.11%. The function of the IM  algorithm about gross error detection came into play when the gross error magnitude was 1.6m. As for gross error detection rate, it exceeded 50% when gross error magnitude was 3.4 m, and it reached 100% when the gross error magnitude was 5 m. Therefore, the gross error detection performance is greatly affected by the magnitude of gross error, and the lower bound of gross error detection depends on the false alarm rate P FA , the ranging accuracy σ , and the geometry of sensors network layout. Thus, improving the ranging accuracy σ , and the geometry of sensors network layout can effectively improve gross error detection performance under guaranteeing the favorable P FA .
The Mean Absolute Error (MAE) and Mean Square Error (MSE) were employed as the indexes to evaluate positioning accuracy. MAE and MSE can be formulated as follows: where ( x, y) is the estimated coordinate of target node, the (x,ỹ) is the true coordinate of target node, N is the total number of epoch. Fig.10 shows the correlation that the IM eliminated the detected single gross error and improved positioning accuracy. According to Fig.10 that when the gross error magnitude exceeding 3 meters, IM made a significant improvement of original positioning accuracy, and as the gross error magnitude increased, the amount of positioning accuracy gradually increased. The MAE improvement was 3.044 m, and the MSE improvement was 12.06 m 2 when the gross error magnitude was 10 m.
Next, we analyzed the performance of IM for two gross error detection in a simulation experiment. Its simulation condition is the same as the single gross error detection. The values of two gross errors were set equal during the simulated experiment process, and the probability of gross errors in each observation was equal. The simulation results are shown in Fig.11 and Fig.12.   FIGURE 11. Comparison between gross error magnitude and two gross error detection rate and identification rate.
According to Fig.11 that when the gross error magnitude was 2.6 m, IM started to detect gross errors, the detection rate reached 50% when gross error was 4.2 m, and the gross error detection rate reached 100 % when the gross error was 8 m. Compared with IM detection of single gross error, the lower bound of two gross errors detection was increased to some extent. In terms of gross identification rate, the minimum identification rate of the IM was 92.67%, the average identification rate of the IM was 98.21%, and the identification performance was slightly reduced relative to single gross error. Fig.12 shows that IM eliminated the detected two gross errors and improved positioning accuracy.
By analyzing the positioning accuracy improvement of IM as in Fig.12, it is found that when gross error was 2 to 4.6 m, the positioning accuracy was reduced since geometry of sensors network deteriorated. The maximum MAE reduction was 0.26576 m and the maximum MSE reduction was

V. EXPERIMENTAL EVALUATION
In this section, a series of field experiments included groups of single point positioning and continuous positioning tests were performed to evaluate the performance of the algorithm proposed.

A. EXPERIMENTAL SETUP
The experimental scene was an underground parking lot of the information department of Wuhan University, whose total area is 2689 m 2 . The underground parking lot as shown in Fig.13. The red rectangles mark the BLE beacons in Fig.13. We can see many bearing columns in this Figure, its height is 2.8m. The intervals among the BLE beacons were properly set in terms of the construction layout and road status of the parking lot. The BLE beacons were installed in the vertical direction on the top of the parking lot at the same height (approximately 3 m).
To clearly illustrate the experiment scene, a digital map of experimental scene was made by the author, as shown in Fig.14.
The red stars represent the BLE beacons location in Fig.14, the total number of BLE beacons was 65. We developed an Android application for collecting positioning data and analyzing data with MATLAB (R2019a, The Math Works, Nitick, MA, USA). The device involved in the experiments was a smartphone running on the Android 8.   beacons can be deployed evenly and non-evenly. According to the simulation analysis in Section IV, the HDOP average value will be minimized when the BLE beacons are evenly deployed, and the positioning accuracy is higher in the same range accuracy. Therefore, an evenly BLE beacon deployment was selected in this work.
The advantages of this deployment are twofold: first, indoor positioning generally only considers the twodimension plane coordinates, the distance in the vertical direction can be seen a constant when the indoor height was known; secondly, the BLE beacons signals propagation is extremely susceptible to obstacles, the ranging value of the receiving end is usually affected by NLOS errors when moving objects and fixed facilities are in the propagation direction, resulting in a decrease in positioning accuracy. The BLE beacons are placed on the top of the room can to ensure the maximum signal propagation in the line-of-sight direction.

B. RSSI RANGING MODEL FITTING
The RSSI collection experimental scenario is shown in Fig.15. The vertical distance between the smartphone end and the Bluetooth node was maintained about 1.805 meters. In the experiment, the horizontal distance between the smartphone and the BLE beacon was varying, and 20 sampling points were set at 0.5-meter intervals from 0 to 9.5 meters. Each sampling point collected 50 sets of RSSI values.  The variation of RSSI at a single sampling point was analyzed. Taking the first sampling point with a horizontal distance of 0 m as an example, its RSSI change curve is shown in Fig.16. As can be seen from Fig.16, between 1∼50 sampling epochs, RSSI is fluctuating between −49db ∼ −57db, with no significant rules. The BLE beacon will send multiple RSSI data within a certain sampling time interval, so it is necessary to average the RSSI data to obtain the ranging value. Therefore, the mean RSSI of each sampling point should be counted while fitting one curve, ranging model so that the model parameters can be calculated. The statistical results of the mean RSSI of the sampling points are shown in table 3. According to Table 3 that as the slant distance increased, the average RSSI decreased overall. Using the logarithmic model in Section III.A for distance-RSSI curve fitting, the logarithmic attenuation model can simplify as shown the formula (36) when the reference distance d 0 is equal to 1 m.
The sampling average results and the data fitting results are visualized in Fig.17. According to Fig.17 that with an increase in the measurement distance, the RSSI value generally shows a downward trend. Between the experimental distance of 1.805 meters and 9.700 meters, the Bluetooth signal strength was attenuated by a total of 26.1 dB, and the attenuation trend conformed to a logarithmic model. After calculation, the unknown parameters a = −31.1811, b = −50.0608, in the formula (36).  Fig.18.
In our research, a two-dimensional plane Cartesian coordinate system was selected. Its origin is the intersection of BLE beacons (outermost) deployed from south to north and the bottom of the parking space (outermost) west to east in Fig.18. The south to north direction is the y-axis and the west to east direction is the x-axis. The real coordinate of the target point A was (19.3689, 27.2417), indicated in Fig. 18, by the black triangle. The BLE beacons were deployed on both sides of the road systematically along the x-axis. Vehicles and pillars have little influence on the Bluetooth signal occlusion. The simulation experiment in Section IV has shown that the number of BLE beacons, n, will affect the IM performance. Therefore, in the field experiment, the number of BLE beacons, n, involved in calculation were set to 8, 10, 12, and 100 sets of positioning data were collected for each n. Scatter diagrams of the positioning results are shown in Fig.19. Black dots denote the positioning results solved only by LSBT without performing IM; Green dots are the positioning results with IM, after MLEGET based on the LSBT, a with recalculated positioning coordinates. It can be seen from the scatter plot of positioning results that the positioning results after performing the IM are densely distributed around the true value point, while the distribution of positioning results before IM was relatively discrete. As the number of BLE beacons, n, increased, the coincident points of the positioning results before and after performing the IM became fewer, that is, the number of increasingly detected gross errors. In addition, the initial positioning results were more discrete, and the improvement of the positioning accuracy of the IM was more apparent.
The absolute positioning error was calculated according to the following formula: where ( x, y) is the estimated coordinate of the target node, the (x,ỹ) is the true coordinate of the target node. The change curve of absolute positioning error before and after the IM is shown in Fig.20. It is can be seen from the variate curve of the absolute positioning error in Fig.20 that the positioning error range fluctuates widely before implementing the IM. The range of fluctuation was more than 5 m, which manifests extreme instability in the signal from the BLE beacons. The fluctuations in the absolute positioning error were smaller after performing the IM, within a range of 3 m. The absolute positioning error curve after performing the IM was below the absolute positioning error curve before the IM, indicating the IM effectively improved positioning accuracy. The cumulative distribution function (CDFs) of the test point A in the case of a different number of BLE beacons are plotted in Fig.21.
According to Fig.21, the positioning accuracy of the LSBT +IM is higher than only LSBT used. The 90% errors of the LSBT-A8+IM, LSBT-A10+IM, and LSBT-A12+IM were 1.9m, 1.8m, and 1.4m, the 90% errors of the LSBT-A8, LSBT-A10 and LSBT-A12 were 2m, 4.1m, and 4.6m, which are reduced by 5%, 56.1%, 60.87% after employing the IM proposed. Also, from Fig.21 it is seen that when the number of BLE beacons is more, the higher the IM performance.
The statistical results using the same evaluation criterion as in Section IV simulation, for test point A are shown in Table 4 and Fig.22. These were MAE, the MSE, and the gross error detection rate(GEDR). Table 4 shows the GEDR by the number of beacons.
It can be seen from Table 4 that as the number n of BLE beacons increased, the GEDR gradually increased (where the GEDR indicated the percentage of gross errors detected    in 100 data sets). The GEDR reached 95% when n = 12.
The RSSI values were sorted in descending order first, and the nearest BLE beacons were selected to participate in calculation. Therefore, the gradual increase in the GEDR also reflects that larger range value, hence the BLE signal is more unstable.
According to the Fig.22, the LBST + IM positioning result is optimal when n = 12, where MAE was 0.90309 m, MAE improvement was 2.280 m, MSE is 0.991 m2, and MSE improvement was 11.354 m2.
The availability of the IM will determine whether to alert the users. Thus, the IM availability was analyzed further to verify the alert mechanism of the IM algorithm in case positioning accuracy no longer meet users' requirement. The HPL represents the maximum positioning error that may be caused by a single gross error that does not reach the specific missed detection rate during the positioning process. Here, the statistical HPL data is compared with the absolute positioning error after performing IM, and the results are plotted in Fig.23, and the statistical results are shown in Table 5.
It can be seen from Fig.23 that when the n is 8, 10, and 12 respectively, the HPL curve can be kept above the positioning error curve. Therefore, the HPL can play a role in protecting the positioning error of BLE RSSI ranging. Comparing HPL to a horizontal alarm level (HAL) can determine whether an alarm is required. Take HAL = 4 m as an example, if the HPL ≥ HAL = 4, the user will be alerted in time.
When the number n of BLE beacons increased from 8 to 12, then the average HPL decreased from 2.8128 to 1.9453; thus, the maximum positioning error that can be generated by the gross error which does not reach the missed detection rate will decrease as the number n of BLE beacons increased.

2) Single point positioning IM analysis in case of different location points
The locations of measure points B, C, D, and E are shown in Fig. 24.
Since the BLE beacons are evenly deployed along the parallel lines on both sides of the road, each section of the road can be regarded as a positioning scene. As the VOLUME 8, 2020     The intervals are 1 meter among these the points in the y-axis direction, as the west side of the E is close to the load-bearing column. When collecting positioning data, the number of BLE beacons was set to 8. At each of these points 100 sets of data were collected. The distribution of positioning results is shown in Fig.25.
According to Fig.25, we analyzed the improvement effect of the IM on the positioning results. The point B near the center of the road displayed large deviations with the most pronounced improvement of positioning results, since the geometry of BLE beacons network was optimal in point B. The initial positioning results of C and D points are close to the positioning results after performing the IM. The positioning results of the E point located at the edge of the positioning scene were severely blocked and discretized to some extent. The statistical results of the points to be measured and the absolute positioning error change curve are shown in Fig.26 and 27.
According to the results displayed in Fig.26, the positioning accuracy was improved after performing the IM to some  extent, in case of points B, C, and D. The positioning accuracy slightly decreased in the point E due to the poor geometry of BLE beacons networks. The overall positioning accuracy has been improved substantially after performing the IM in these test points experiments as can be seen from Fig.27.
The CDFs of positioning error of test points B, C, D, E are plotted in the figure28. According to Fig.28 that when only the LSBT method was performed, the maximum value of positioning results was close to 6m, but when LSBT+IM was performed, the maximum value of the positioning result not exceeding 4m. The error at a 95% confidence level LSBT+IM-B, was 2.2m, for LSBT+IM-C the error was 1.8m, for LSBT+IM-D the error was 2.0m, and for LSBT+IM-E the error was 2.3m. The errors at a 95% confidence level for LSBT-B was 2.8m and 2.5m for LSBT-C. For LSBT-D the error was2.5m, and for LSBT-E the error was 1.9m, which were reduced by 21.43%, 28.0%, 20%, and -21.5% after performing the IM proposed.
The HPL of measured points and the absolute positioning error after performing the IM are shown in Fig.29. It can be seen from the HPL changing graph that changing the position of measure points and the degree of occlusion, very few positioning errors are greater than HPL due to random fluctuations in signal strength or gross error mistake identification. HPL can still guarantee positioning in most of the time above the error line, it illustrates the HPL can protect the positioning results.

2) CONTINUOUS POSITIONING RESULTS ANALYSIS
To further verify the performance of the proposed IM during the continuous positioning process, we conducted 4 groups continuous positioning test experiments. The true path in continuous positioning is depicted in Fig.30 during the experiment. The total length of the true path is 42 m, the lengths across y-axis and x-axis directions are 19 m and 23 m respectively.
The experimenter held the smartphone to capture a certain degree of signal occlusion so that the observations containing gross errors. A speed of movement was constant during the experiment at about 0.75 m per second. During the experiment, four test datasets were collected in the forward and reverse directions along the true path. Each dataset contained 45 samples. The trajectories of the collected data points are shown as a plot in Fig. 31.
The proposed IM algorithm significantly improved continuous positioning accuracy. It can be seen from Fig.31 that four sets of LBST test experiments have large positioning errors due to the effect of RSSI value fluctuations. The points that deviate from the true path mostly appeared while traveling along the x-axis direction. In test experiment A, IM optimized the LBST trajectory slightly. In test experiments, B, C, and D, discrete points that deviate from the true path were effectively corrected after performing IM. The positioning points after performing the IM were within 3 meters of the true path, thus test route was closer to the real route, and the maximum of the LBST was 6m.

VI. DISCUSSION AND CONCLUSIONS
Innovatively, this article has proposed an integrity monitoring (IM) algorithm based on the parity vector and Maximum Likelihood Estimation Gross Error test (MLEGE) for indoor positioning using Bluetooth Low Energy (BLE) beacons. The proposed IM algorithm has achieved significant improvements on the BLE positioning accuracy and positioning robustness. This IM algorithm can identify and eliminate two gross errors simultaneously to improve positioning accuracy during the positioning process; most important of thing, the IM provides one technique of timely alert when the positioning system is unavailable due to the positioning accuracy no long satisfies user requirement. A series of simulation and field experiments were conducted to evaluate the performance of the IM algorithm. The discussion and conclusions about simulation analysis and field experiments as follows.

A. SIMULATION RESULTS AND ANALYSIS
According to the simulation analysis for the impact of Horizontal Dilution of Precision (HDOP) on Horizontal Protection Level (HPL), the impact of HDOP on gross errors detection performance and the impact of the different magnitudes of gross errors detection performance in the Section IV, we can conclude that the gross errors detection performance of proposed IM algorithm mainly affected by geometry of BLE network, magnitude of gross errors, false alarm rate and ranging accuracy.
When the BLE beacons were evenly deployed, the HPL gradually decreased as the number of BLE beacons increased, which meant the positioning accuracy will be improved in case of the number of BLE beacons increased during the positioning process.
In 1000 times of simulation experiment, the average value of gross error detection rate for single gross error and the two gross errors were 98.9% and 85.9% respectively in case of the number of BLE beacons was 6 and magnitude of the gross error was 5 meters. It's illustrated that the IM can effectively detect and identify gross errors. According to the analysis for the impact of HDOP on gross errors detection performance can be concluded that as the HDOP increased, the gross error detection rate gradually decreased, and the identification rate of gross error gradually increased.
According to analysis for the Section IV.C, the gross error detection performance is greatly affected by the magnitude of gross error, and the detection performance on single gross error was slightly higher than that on two gross errors. For single gross error, the IM algorithm gross error detection came into play when the gross error magnitude was 1.6m. The gross error detection rate reached 100% when the gross error magnitude was 5 m; For two gross errors, IM started to detect gross errors when gross errors magnitude was 2.6 meters, and the gross error detection rate reached 100 % when the gross error magnitude was 8.0 meters.
When the gross error exceeding 3 m, the proposed IM can effectively improve the positioning accuracy, and as the gross error magnitude increased, the amount of positioning accuracy improvement also gradually increased.
In a short, the proposed IM algorithm for gross errors detection is working effectively.

B. FIELD EXPERIMENT RESULT AND ANALYSIS
Point positioning and continuous positioning field experiments were conducted to verify the performance for the proposed IM, Field experiments showed that the IM algorithm supports the positioning accuracies of <1.9143 m at 90% of time (average of 7 groups single points), which performed better than <2.9143 m of the least square based on Taylor expansion (LSBT); For the continuous positioning (4 groups test experiments), the maximum error of the IM not exceeded 3 m, which performed better than 6 m of the LSBT method.
In the single points positioning experiments with a different number of BLE beacons (i.e. in test point A). The 90% errors of the LSBT-A8+IM, LSBT-A10+IM, and LSBT-A12+IM were 1.9m, 1.8m, and 1.4m, the 90% errors of the LSBT-A8, LSBT-A10 and LSBT-A12 were 2m, 4.1m, and 4.6m, which are reduced by 5%, 56.1%, 60.87% after employing the IM proposed. Therefore, we can conclude that as the number of BLE beacons is more, the higher the IM performance.
According to the single points positioning experiments in case of different location points (i.e. in test points B, C, D, E), we can conclude that the positioning accuracy after performing the IM was most higher in test point B while the positioning accuracy after performing the IM was most lower in the test point E. Since the geometry of BLE beacons network was optimal in point B which near the center of the road. By contrast, the geometry of BLE beacons networks was poor in point E due to the effect of the load-bearing column in the parking lot on the Bluetooth signal occlusion became stronger. Thus, the geometry of BLE beacons network would be affected the performance of the IM, the geometry of BLE beacons network more optimal, the higher the IM performance.
In case of positioning accuracy, no longer matches users' expectation, the alert technique of the IM algorithm was further verified. The experimental results showed that the users will be alerted in time to avoid risk from positioning accuracy no longer satisfying the user demand. This is an indispensable function of the positioning system, particularly in firefighting services, peacekeeping missions, emergency safety, disaster relief, mobile health, and intelligent robot applications. where n represents the aggregate of BLE beacons observations participating in the positioning at the current time, and u is the necessary number of BLE beacons observations for calculating the coordinate estimation value, Q is an identity orthogonal matrix of order n × n, and R is (n × u) an upper triangular matrix. Bring (A-2) into equation (A-1) and left multiply both sides by Q T can be obtained: where U is an (n × n) identity orthogonal matrix, the is a (n × u) matrix, V is an (u × u) identity orthogonal matrix.