Application of Baro-Altimeter Sensor in Emergency Positioning and Navigation Based on Compass GEO Satellites

Abstract

Frame synchronization time of geostationary Earth orbit (GEO) satellites is faster than that of incline geosynchronous satellite orbit and medium Earth orbit satellites in Beidou navigation satellite system (BDS). So time to first fix can be reduced greatly with GEO satellites (GEOs) to provide emergency positioning service during the cold start process. Meanwhile, GEOs can also be utilized to provide radio determination satellite services. Therefore, it is meaningful to research the emergency positioning algorithm based on GEOs. However, due to the poor geometry of the GEOs, the availability and accuracy degrade distinctly. To solve this problem, a low-cost high resolution MS5534B barometric (baro-) module is applied to BDS. Firstly, the principle of emergency positioning based on a baro-altimeter sensor and its performance such as the accuracy are elaborated. Then the effects of baro-altimeter sensor measurement error on positioning are analyzed. Finally, after analyzing the limitation of the conventional algorithms, a new high-accuracy emergency positioning algorithm with baro-altimeter sensor aiding is proposed, which is not limited by the integration and user's altitude. In comparison to the conventional algorithms, the stimulant and experimental results validate the effectiveness, robustness, and higher performance of the proposed algorithms incorporated with sensors.

1. Introduction

Beidou-2 regional satellite navigation system (BDS) has been providing accurate and continuous three-dimensional (3-D) position, communication, and time services since 27th, December, 2012. However, in the typical urban environment, the signal blockage caused by the tall buildings reduces the satellite availability and deteriorates the observation geometry, with the extreme case being solution unavailability [1]. Even if the number of visible satellites is more than four, the geometry is poor; it is still unable to obtain the correct positioning result. To overcome the system weakness, a method that combines BDS with sensors or pseudolites can be adopted to provide emergency positioning service [2]. For BDS, the emergency positioning is of great significance for increasing the availability of radio navigation satellite system (RNSS) and radio determination satellite system (RDSS).

Among the three types of satellite in RNSS navigation subsystem, the geostationary Earth orbit (GEO) satellite has two conspicuous features compared to inclined geosynchronous satellite orbit (IGSO) and medium Earth orbit (MEO) satellites. The first is that the coverage region is fixed and the satellites are always visible in most of the Asian Pacific region. The second is that the broadcasted D2 navigation message rate is 10 times of the rate of D2 navigation message broadcasted by MEO/IGSO, and its rate can reach up to 500 bps [3]. That is, in the same condition, the frame synchronization time of the navigation message is 10.8  seconds shorter than that of MEO/IGSO. Therefore, for the user in Asian Pacific region, when in the cold start situation, the acquisition strategy which takes GEO satellite as the first acquisition satellites can reduce the time to first fix time (TTFF) by about 10 seconds. Here, the precondition is that the emergency positioning must be used so as to overcome the influences of the poor geometry of the GEO constellation on the positioning accuracy result.

RDSS navigation subsystem is the inheritance and development of the Beidou-1 positioning system, and it is an active positioning system consisting of five GEO satellites which are shared by RNSS. Nonetheless, due to the disadvantages of poor disguising, limited service area, and large error geometric coefficient, the application scope of the system are strongly limited. To overcome the defect of the active positioning system, many kinds of emergency positioning methods have been proposed in recent years. In literatures [4, 5], a positioning method is proposed by combining the Beidou three satellites with other navigation system. Compared to a single system, the dilution of precision (DOP) and positioning accuracy are both improved. In literature [6], a new scheme is proposed by applying two pseudolites to augment RDSS, and an optimum distribution in the augmentation condition with two pseudolites is given. However, due to the sight distance, its operational domain is limited in order to guarantee the positioning accuracy. In literature [7], a hyperbola positioning system solution was proposed by adding intersatellite signal transmission links to realize intersatellite differential. Since the positioning error is amplified, the modification of the system and the receiver is large. Further, based on the method of GPS baro-altimeter altitude aiding [8], an improved passive positioning algorithm is proposed by combining Beidou three satellites and barometric altitude observation equation [9, 10]. The coverage area of this scheme is wider than that of pseudolites, and it has the advantage to overcome the problem that the visible satellites are not enough due to high mask angle.

Even now, previous research still has some shortcomings as follows. (1) The pseudorange error and altitude error are always not complying with the same distributions, which will influence the reliability of positioning accuracy evaluation. (2) There is some linear error between the initial model and the altitude observation model obtained by Taylor expansion about the approximate position. (3) The emergency positioning algorithm based on the conventional observation models has limitation of being unable to compute the position fix when the user's altitude exceeds over about $3.6e7\hspace{0.17em}\text{m}$ and there are more than 3 GEO satellites in the integration. (4) The lack of GEO satellite number is not the essential reason of unavailability, and the emergency positioning is also needed even if the number of visible GEO satellites is more than three in BDS. However, the literatures above only considers the situation of three satellites in the integration. In terms of the integration of GEO satellites with 4 or 5 GEO satellites, the performance of the conventional aiding algorithm will change, even reduce, which will be discussed below.

It is of significant practical interest, therefore, to enhance the accuracy and availability of GEO-based RDSS and RNSS subsystem using an external source of baro-altimeter sensor's altitude information in emergency environment. The objective of the work presented here is to introduce a novel station-centric-intersection emergency positioning algorithm by incorporating a proposed error-free altitude observation model into the basic GEO measurement equation set and weighting of pseudorange and baro-altimeter measurements, instead of the conventional algorithms. We show that improvements of the accuracy, availability, and robustness can be achieved with the proposed algorithm proposed. This paper is organized as follows. In Section 2, GEO satellites and baro-altimeter sensor integration principle is elaborated, along with the sensor accuracy budget. Section 3 analyzes the modeling error of the three conventional altitude observation models and proposes an accurate altitude observation model. Section 4 introduces how the altitude observation constraint can be incorporated into the basic equation for the set of BDS GEO measurements and demonstrates its contribution to emergency positioning accuracy. Section 5 presents the performance (as measured by coverage characteristic, accuracy, and robustness) of the algorithm proposed with simulation and observation data, in comparison to the conventional algorithm, along with the estimation of GEO user equivalent range error (UERE). Finally, Section 6 presents the conclusions and future work.

2. GEO Satellites and Baro-Altimeter Sensor Integration Principle and Accuracy Budget

2.1. GEO Satellites and Baro-Altimeter Sensor Integration Principle

RDSS navigation subsystem is a transmitting navigation system consisting of ground measuring and control center (MCC), calibration system, GEO satellites, and user receiver. The GEO satellite retransmits the outbound signal generated from MCC to the coverage area, and the user receiver measures the pseudorange by calculating the propagation delay of the signal transmitted from MCC to the receiver to computer the position fix in the passive location mode. Therefore, the receiver must make at least four pseudorange measurements to satellites.

Consider

$$\begin{array}{cc}& {\rho }_i=\sqrt{{\left(x_i-x_u\right)}^{\mathrm{2}}+{\left(y_i-x_u\right)}^{\mathrm{2}}+{\left(z_i-z_u\right)}^{\mathrm{2}}}+c{t}_u+{\epsilon }_i,\end{array}$$

(1)

where ${\rho }_i$ is the pseudorange measurement from the receiver to the ith ( $i\ge 4$ ) satellite; ( $x_i,y_i,z_i$ ) is the unknown 3D position in ECEF coordinates when the pseudorange measurement was taken; c is the speed of light; $t_u$ is the receiver clock offset; and ${\epsilon }_i$ is the pseudorange measurement error.

In the canyon or high-rise building environment, due to the signal blockage or the limitation of beam coverage, the phenomenons that less than four satellites are available and the current geometry of satellites is poor often occur. As a result, the direction cosine matrix is rank deficient or singular, which leads to a failure of computing the position fix. A method using additional information from an external source can be adopted to obtain one or two unknown components of the user's position. In principle, it is equivalent to adding one or two available satellites in the constellation, and therefore, the emergency position fix can be realized. Generally, eliminating the unknown local clock error or the unknown altitude in geodetic coordinate reference frame is the two methods that could be used. However, eliminating the local clock error needs strict time synchronization between the user receiver and BDT, which is very difficult in engineering. In addition, if the 3D unknown position could be reduced to 2D, the problem can be solved. In terms of the 3D position components, the vertical one could be obtained relatively easy, which can be obtained from a baro-altimeter sensor. Since a baro-altimeter sensor is not limited by the region and geographic position, the availability is complemented and improved, and the application scope is wide. Figure 1 illustrates the emergency positioning principle, and (2) gives the measured baro-altimeter altitude equation in BDS.

Figure 1

Emergency positioning principle based on 5 GEO satellites and baro-altimeter sensor in BDS.

Consider

$$\begin{array}{cc}& h=f\left(x_u,y_u,z_u\right)+{\epsilon }_h,\end{array}$$

(2)

where h is the geodetic height in China geodetic coordinate system 2000 (CGCS2000) from the measured baro-altimeter altitude [11], $f(x_u,y_u,z_u)$ is the altitude observation model, ${\epsilon }_h$ is the error in measured baro-altimeter altitude after it has been calibrated. To simplify the analysis, here, h is assumed to have calibrated all the modeling error sources, such as the height difference between the mean sea level and the geoid [12].

2.2. Accuracy Budget of Baro-Altimeter Sensor

Similar to the pseudorange measurement, the baro-altimeter altitude measurement is also corrupted by errors. From user's point of view, what they are concerned with is that in the emergency positioning, which kind of effect the baro-altimeter altitude error will bring to the positioning accuracy. That is, it is necessary to consider whether the specification of the baro-altimeter can meet the system accuracy requirement.

Thus, according to the model provided in [12], the effects of the error ${\epsilon }_h$ on the positioning error is figured out with the integration of three GEO satellites with a baro-altimeter sensor, as shown in Figure 2.

Figure 2

(a) Effects of the standard deviation of the baro-altimeter measurement error on horizontal and vertical positioning errors. (b) Microbaro-altimeter module MS5607-02BA (Intersema Inc.).

Given the user geodetic coordinates ( $\phi ,\lambda ,h$ ), Figure 2(a) shows the variation of the horizontal positioning error and vertical positioning error with the baro-altimeter altitude error ${\epsilon }_h$ , at $\phi =15$ °N, $\phi =35$ °N, and $\phi =55$ °N, respectively, while the longitude λ and the altitude h keep constant, $\lambda =110$ °E, $h=10$  m. It can be seen that the horizontal positioning error is monotonously increasing with ${\sigma }_h$ , and the larger the latitude, the bigger the horizontal positioning error. Nonetheless, the vertical positioning error is proportional to ${\sigma }_h$ , with the slope $k_{\text{slope}}=1$ and is independent on λ. Therefore, in the emergency positioning, the error of the baro-altimeter altitude sensor will directly affect the user positioning accuracy. For example, assuming that the user is located at (35°N, 110°E, and 100 m), to meet the requirement of the horizontal positioning error of 100 m and the vertical positioning error of 10 m, the measured altitude error must be smaller than 54.9 m.

There are various types of baro-altimeter altitude sensor products, and the series products based on MEMS technology are widely used in personal navigation devices, such as the MS5534B and MS5607-02BA03. As the temperature compensation algorithm and filtering algorithm are used, the effects of the temperature variation and barometric turbulence are reduced, and the measurement accuracy is remarkably increased. Further, test result shows that under nominal standard atmospheric conditions, the average standard deviation of relative altitude measurement with MS5534B is 0.2 m [13], which completely fulfils the accuracy requirement in emergency positioning [14]. Figure 2(b) shows a product picture of MS5607-02BA [15].

3. Altitude Observation Model

3.1. Error Analysis and Limitation of Conventional Altitude Observation Model

When less than four pseudorange measurements are available, (1) cannot be solved to accomplish position fix. According to the aided principle, if the baro-altimeter altitude observation equation (2) is added to the basic GEO measurement equation; theoretically, the unknown user position and clock error can be solved. However, the geodetic height h is an independent parameter in geodetic coordinate reference frame, and the unknowns $x_u$ , $y_u$ , and $z_u$ in (1) are three independent parameters in ECEF. Obviously, the two coordinate systems are not orthogonal to each other. Thus, h is dependent on $x_u$ , $y_u$ , and $z_u$ , simultaneously. In other words, the effects of the error of h are not only on vertical accuracy itself, but also are on the horizontal accuracy. Therefore, how to model the measured baro-altimeter altitude accurately is the key point. Conventional observation models have been proposed for the measured baro-altimeter altitude, but due to using approximations for various reasons, the modeling error is induced, which would affect the performance of the algorithm. For simplification, this paper focuses on three models adopted in many papers, which are described as follows.

3.1.1. Approximate Ellipsoid Model

Figure 3(a) depicts the first conventional model called approximate ellipsoid model, which is adopted in many papers. In the figure, a and b are, respectively, the semimajor axis and the semiminor axis of the Earth ellipsoid, while $\mathrm{\Delta }a$ and $\mathrm{\Delta }b$ are, respectively, the semimajor axis variation and the semiminor axis variation of the approximate ellipsoid compared to the reference ellipsoid. Assume that $P_1(x_1,y_1,z_1)$ is one point in reference ellipsoid, and $P_2(x_u,y_u,z_u)$ is the user position point of which the geodetic height is h. The relationship between $P_1$ and $P_2$ is that $P_2$ lies above the reference ellipsoid at an altitude h in the direction of the ellipsoid normal at the user's receiver that passes $P_1$ point. Nonetheless, the approximate ellipsoid model considers the user position point $P_2$ is located on the ellipsoid that meets the condition $\mathrm{\Delta }a=\mathrm{\Delta }b=h$ , which can be expressed as

$$\begin{array}{cc}& \mathrm{1}=\frac{{\left(x_u\right)}^{\mathrm{2}}}{{\left(a+h\right)}^{\mathrm{2}}}+\frac{{\left(y_u\right)}^{\mathrm{2}}}{{\left(a+h\right)}^{\mathrm{2}}}+\frac{{\left(z_u\right)}^{\mathrm{2}}}{{\left(b+h\right)}^{\mathrm{2}}},\end{array}$$

(3)

where h is an unknown parameter; $(x_u,y_u,z_u)$ is an unknown position vector; However, $P_2$ is truly located on the curved surface which is formed by all the points that have the same height h above the reference ellipsoid, denoted as $S_1$ . Therefore, the coordinates solved by (3) is the approximation for $P_2$ . That is to say, the model in (3) includes modeling errors. Furthermore, to evaluate the effect of the error, assume that there is $h^{\prime }$ that meets the ellipsoid equation in (4). Consider the geometric constraint relationship of $P_2$ and $P_1$ , given $(x_u,y_u,z_u)$ , and the following equations can be obtained:

$$\begin{array}{cc}& \mathrm{1}=\frac{x_u^{\mathrm{2}}}{{\left(a+h^{\prime }\right)}^{\mathrm{2}}}+\frac{y_u^{\mathrm{2}}}{{\left(a+h^{\prime }\right)}^{\mathrm{2}}}+\frac{z_u^{\mathrm{2}}}{{\left(b+h^{\prime }\right)}^{\mathrm{2}}},\end{array}$$

(4)

$$\begin{array}{cc}& {\mathbf{X}}_u=\left[\begin{array}{c}{x}_u\\ y_u\\ z_u\end{array}\right]=\mathbf{F}{\mathbf{X}}_{\mathrm{1}}=\left[\begin{array}{ccc}\mathrm{1}+\frac{h}{u{a}^{\mathrm{2}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}+\frac{h}{u{a}^{\mathrm{2}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}+\frac{h}{u{b}^{\mathrm{2}}}\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{1}}\\ y_{\mathrm{1}}\\ z_{\mathrm{1}}\end{array}\right],\end{array}$$

(5)

$$\begin{array}{cc}& \mathrm{\Delta }h=h^{\prime }-h=f_e\left(x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}}\right),\end{array}$$

(6)

where $h^{\prime }$ is an unknown parameter, and $(x_u,y_u,z_u)$ is a known position vector. ${\mathbf{X}}_{\mathrm{1}}$ and ${\mathbf{X}}_u$ in (5) are the column vectors of $P_1$ and $P_2$ , respectively. F depicts the relationship between $P_1$ and $P_2$ , which is a diagonal matrix, and μ is the norm of vector ( $x_1/a^2,y_1/a^2,z_1/a^2$ ). Given ( $x_u,y_u,z_u$ ) in (4), by combining (4) and (5), the characteristic of modeling error $\mathrm{\Delta }h$ can be obtained, as shown in Figure 4.

Figure 3

(a) Approximate ellipsoidal model. (b) Approximate ellipsoid model errors at different user's position coordinates ( $x_u,y_u,z_u$ ). (c) Approximate ellipsoid model modeling errors at different h.

Figure 4

(a) Earth-centered fictitious satellite model. (b) θ of the Earth-centered fictitious satellite model at different latitude φ and longitude λ when $h=100$  m.

Figure 3(b) depicts $\mathrm{\Delta }h$ at different geodetic height h. To obtain the maximum magnitude of the error $\mathrm{\Delta }h$ , θ is defined, to go through all position points of $P_1$ on the reference ellipsoid, by using the parametrizations $x^2+y^2=a^2{\cos }^2(\theta )$ and $z^2=b^2{\sin }^2(\theta )$ . Figure 3(c) shows the variation of the modeling error $\mathrm{\Delta }h$ with geodetic height h. Here, at each h, $\mathrm{\Delta }h$ is the maximum modeling error value over all the θ. It can be seen that $\mathrm{\Delta }h$ is proportional to h. Thus, when the geodetic height is 2000 m, the maximum induced modeling error is 0.028 m, which cannot be ignored in high-accuracy positioning applications.

3.1.2. Earth-Centered Fictitious Satellite Model

Figure 4(a) depicts the second conventional model called Earth-centered fictitious satellite model [12, 16]. The model views the baro-altimeter altitude as a fictitious satellite at the center of the Earth, which, in principle, is equivalent to adding one pseudorange measurement between the user and the fictitious satellite. Equation (7) is the corresponding altitude observation equation. According to the altitude definition in CGCS2000, the geodetic height is strictly defined in the ellipsoid normal at the user's receiver passing $P_2$ . However, since the reference ellipsoid is not a sphere, obviously, the angle w between the ellipsoid normal and the line-of-sight pseudorange is not zero. It has been shown that when the receiver position is (45°N, 0°E, and 100 m), within the range of $\pm 100$  m at the fictitious sphere surface, the maximum distance error between the fictitious sphere and the reference ellipsoid is 0.4 m [12]. Furthermore, this error is approximately proportional to $\sin w$ , and w is dependent only upon the receiver position, which is shown in Figure 5.

Figure 5

Normal fictitious model.

Consider

$$\begin{array}{cc}& {\rho }_v=\sqrt{{\left(x_v-x_u\right)}^{\mathrm{2}}+{\left(y_v-y_u\right)}^{\mathrm{2}}+{\left(z_v-z_u\right)}^{\mathrm{2}}}+{\epsilon }_h,\end{array}$$

(7)

where ${\rho }_v$ is the fictitious pseudorange between the fictitious satellite and the user after the altitude is converted, ( $x_v,y_v,z_v$ ) is the 3D position coordinates of the fictitious satellite.

It can be seen from Figure 4(b) that w is dependent on φ, which is almost not affected by λ. Within the system coverage range (φ: 5~55°, λ: 75°~−155°), when $\phi ={45}^{°}$ , w reaches its maximum, $w_{\max }={0.1924}^{°}$ .

3.1.3. Normal Fictitious Satellite Model

Figure 5 depicts the third conventional model called normal fictitious satellite model. In the model, the fictitious satellite is located at the intersection of the ellipsoid normal at the user's receiver line and z-axis, of which the 3D coordinate is ( $0,\hspace{0.17em}0,\hspace{0.17em}\text{and}\hspace{0.17em}-N{e}^2\sin \phi$ ). Thus, the following altitude observation equation is required:

$$\begin{array}{cc}& h=\sqrt{x_u^{\mathrm{2}}+y_u^{\mathrm{2}}+{\left(z_u+N{e}^{\mathrm{2}}\sin \phi \right)}^{\mathrm{2}}}-N+{\epsilon }_h,\end{array}$$

(8)

where N is the prime vertical radius.

In terms of the model definition, theoretically, the normal fictitious satellite model completely conforms to the definition of the reference ellipsoid coordinates. Therefore, the model itself is error-free. However, when it is linearly expanded in a Taylor series at the approximate position ( $x_{u0},y_{u0},z_{u0}$ ), the unknown fictitious satellite position ( $0,\hspace{0.17em}0,\hspace{0.17em}-N{e}^2\sin \phi$ ) used in the direction cosine matrix is dependent on the unknown $x_u$ , $y_u$ , $z_u$ . Strictly, for Taylor expansion, the parameters N and φ cannot be treated as constants in the linearization process. Thus, the accurate expansion formula is defined as follows:

$$\begin{array}{l}\mathrm{\Delta }h=\frac{x_u+\left(\left(\partial f_z\left(N,\phi \right)\right)/\partial x\right)}{r_h}\mathrm{\Delta }x\\ +\frac{y_u+\left(\left(\partial f_z\left(N,\phi \right)\right)/\partial y\right)}{r_h}\mathrm{\Delta }y+\frac{\left(\partial f_z\left(N,\phi \right)\right)/\partial z}{r_h}\mathrm{\Delta }z,\end{array}$$

(9)

where

$$\begin{array}{cc}& f_z\left(N,\phi \right)=z+N{e}^{\mathrm{2}}\sin \phi ,\\ r_h=\sqrt{x_u^{\mathrm{2}}+y_u^{\mathrm{2}}+{\left(z_u+N{e}^{\mathrm{2}}\sin \phi \right)}^{\mathrm{2}}}.\end{array}$$

(10)

It can be seen from (10) that a closed-form Taylor expansion is hard to achieve. Therefore, to reduce the computational burden of computing the position solution, the approximations are usually adopted in the linearization, as shown in (10) [17]. Obviously, some linear modeling error is induced, which may impact the final positioning result.

Consider

$$\begin{array}{cc}& \frac{\partial f_z\left(N,\phi \right)}{\partial x}=\frac{\partial f_z\left(N,\phi \right)}{\partial y}=\frac{\partial f_z\left(N,\phi \right)}{\partial z}=\mathrm{0}.\end{array}$$

(11)

3.2. Improved Altitude Observation Model

Through the analysis of the conventional observation modeling error, it can be seen that the three models are all modeled in ECEF coordinates. However, the geodetic altitude does not directly correspond to any component of x, y, z in ECEF coordinates. Thus, when we try to establish an error-free model by choosing $\mathrm{\Delta }\mathit{u}=(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)$ in ECEF as the observation vector, corresponding to the altitude observation equation in the direction cosine matrix H will be extremely complex, which may even result in no closed-form expression.

Therefore, to generate an error-free altitude observation model to accurately evaluate the positioning performance, the key point is to choose a proper coordinate system. Considering the user geodetic height is simply the minimum distance between the user and the reference ellipsoid, which is an independent parameter in local tangent plane coordinates (east, north, and up), this paper proposes an improved altitude observation model suitable for RDSS and RNSS passive positioning. In the improved mode, the parameter definition and Newton iterative calculation process are both realized in ENU coordinates. Since the ENU coordinates define the user position point as the origin of the coordinate, in comparison to conventional model, the observation equation and the altitude observation vector by linear Taylor expansion have more simple and intuitive format, as shown in (12).

Consider

$$\begin{array}{cc}& h=z_u+{\epsilon }_h,\end{array}$$

(12)

where subscript IM indicates the observation model used, and $Z_{\text{IM}}$ is the component of the receiver 3D position vector at the z-axis in ENU coordinates. It should be noted that as the altitude observation does not involve with time information, if the receiver displacement in the ENU coordinates is chosen as the state vector of the positioning equation, then the additional altitude observation vector of the direction cosine matrix H can be expressed as

$$\begin{array}{cc}& {\mathbf{H}}_{4,:}=\left[\begin{array}{cccc}0& 0& 1& 0\end{array}\right].\end{array}$$

(13)

It can be seen that in the improved altitude observation model, only $z_u$ component is dependent on the state vector. Even more importantly, the corresponding parameters in the row vector ${\mathbf{H}}_{4,:}$ of the direction cosine matrix H is constant; therefore, the first order Taylor expansion is strict, avoiding the linear error and modeling error. Thus, it can be expected that performance would be enhanced by applying the error-free proposed improved power observation model, compared to conventional models, as will be discussed next.

4. Emergency Positioning Algorithm with Altitude Observation and Contribution

4.1. Altitude-Aided Station-Centric-Intersection Emergency Positioning Algorithm

The positioning accuracy is determined by range errors and user-satellite geometry. Dilution of precision (DOP) reflects the user-satellite geometry and represents the amplification of the standard deviation of the measurement errors onto the solution. Smaller DOP often represents good satellite geometry. Therefore, DOP is one of the most important measures to evaluate the satellite navigation system performance.

Generally, when evaluating the positioning accuracy with DOP, it is always assumed that range errors are zero mean with a Gaussian distribution, and that range errors are independent from satellite to satellite [18, 19]. But this assumption is not always effective. In the passive positioning system with baro-altimeter aiding, the statistics property of the GEO satellite pseudorange and the statistics property of the altitude measurements are not completely the same, which is reflected in the inequality of the satellite standard deviation ${\sigma }_{\text{geo}}$ and the altitude measurement error standard deviation ${\sigma }_h$ . If still evaluating positioning accuracy under the assumption in [7, 10, 17], the obtained position solution will not be the optimum. To solve the evaluation problem of unequal measurement variance, weighted dilution of precision (WDOP) is introduced to replace the conventional DOP. Weighting factor is used to weight each pseudorange measurement according to their relative noise level and the relative importance of each estimation. Thus, the positioning accuracy is more precisely evaluated [16, 20]. In general, this paper refers to a widely used WDOP definition [16, 21], as shown in (14).

Consider

$$\begin{array}{cc}& \text{WDOP}=\sqrt{\mathrm{tr}{\left({\mathbf{H}}^T{\mathbf{W}}^{-\mathrm{1}}\mathbf{H}\right)}^{-\mathrm{1}}},\end{array}$$

(14)

where H is the direction cosine matrix, W is the weighting matrix, and $w_i\hspace{0.17em}\hspace{0.17em}(i=1,\dots ,n)$ is the weighting factors of the ith range error. Assuming that the pseudorange measurement error is ${\sigma }_i^2$ , and has an independent Gaussian distribution, ${\sigma }_{\text{UERE}}$ is the pseudorange error factor. We define weighting factor as shown in (15), which is the ratio of ${\sigma }_i^2/{\sigma }_{\text{UERE}}^2$ .

Consider

$$\begin{array}{l}\mathbf{W}=\mathrm{diag}\left[\begin{array}{llll}{w}_{\mathrm{1}}& w_{\mathrm{2}}& \cdots & w_n\end{array}\right]\\ =\mathrm{diag}\left[\begin{array}{llll}\frac{{\sigma }_{\mathrm{1}}^{\mathrm{2}}}{{\sigma }_{\text{UERE}}^{\mathrm{2}}}& \frac{{\sigma }_{\mathrm{2}}^{\mathrm{2}}}{{\sigma }_{\text{UERE}}^{\mathrm{2}}}& \cdots & \frac{{\sigma }_n^{\mathrm{2}}}{{\sigma }_{\text{UERE}}^{\mathrm{2}}}\end{array}\right]\end{array}$$

(15)

when computing the position fix using a WLS algorithm based on WDOP, since the pseudorange with the larger error contributes less to the positioning solution, the corresponding weighted DOP can reflect the positioning error more precisely. Furthermore, it is noted that DOP can be categorized into the concept system of WDOP, which is the special form of WDOP when $\mathit{W}=\mathrm{diag}\left[\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right]$ .

Taking RDSS 2-D passive positioning with baro-altimeter altitude sensor aiding as an example, according to the analysis above, the set of observation equations in ENU coordinates can be obtained as follows:

$$\begin{array}{cc}& {\rho }_{\mathrm{1}}=\sqrt{x_{\mathrm{1}E}^{\mathrm{2}}+y_{\mathrm{1}N}^{\mathrm{2}}+z_{\mathrm{1}U}^{\mathrm{2}}}+c{t}_u+{\epsilon }_i,\\ {\rho }_{\mathrm{2}}=\sqrt{x_{\mathrm{2}E}^{\mathrm{2}}+y_{\mathrm{2}N}^{\mathrm{2}}+z_{\mathrm{2}U}^{\mathrm{2}}}+c{t}_u+{\epsilon }_i,\\ {\rho }_{\mathrm{3}}=\sqrt{x_{\mathrm{3}E}^{\mathrm{2}}+y_{\mathrm{3}N}^{\mathrm{2}}+z_{\mathrm{3}U}^{\mathrm{2}}}+c{t}_u+{\epsilon }_i,\\ h=z_{hU}+{\epsilon }_h,\end{array}$$

(16)

where ${\rho }_i$ is the pseudorange measurement between a user and the ith GEO satellite; ( $x_{iE},y_{iN},z_{iU}$ ) is the 3D position coordinate in ENU coordinates, and the definition of $h_i$ and $z_{hU}$ is the same as that in (12). After first-order Taylor linear expansion, a direction cosine matrix H can be obtained as follows:

$$\begin{array}{l}\mathbf{H}=\left[\begin{array}{llll}\frac{x_{\mathrm{1}E}}{r_{\mathrm{1}}}& \frac{y_{\mathrm{1}N}}{r}& \frac{z_{\mathrm{1}U}}{r_{\mathrm{1}}}& 1\\ \frac{x_{\mathrm{2}E}}{r_{\mathrm{2}}}& \frac{y_{\mathrm{2}N}}{r_{\mathrm{2}}}& \frac{z_{\mathrm{2}U}}{r_{\mathrm{2}}}& 1\\ \frac{x_{\mathrm{3}E}}{r_{\mathrm{3}}}& \frac{y_{\mathrm{3}N}}{r_{\mathrm{3}}}& \frac{z_{\mathrm{3}U}}{r_{\mathrm{3}}}& 1\hspace{0.17em}\\ 0& 0& 1& 0\end{array}\right]\\ =\left[\begin{array}{llll}-\cos {\theta }_{\mathrm{1}}\sin {\alpha }_{\mathrm{1}}& -\cos {\theta }_{\mathrm{1}}\cos {\alpha }_{\mathrm{1}}& -\sin {\theta }_{\mathrm{1}}& \mathrm{1}\\ -\cos {\theta }_{\mathrm{2}}\sin {\alpha }_{\mathrm{2}}& -\cos {\theta }_{\mathrm{2}}\cos {\alpha }_{\mathrm{2}}& -\sin {\theta }_{\mathrm{2}}& 1\\ -\cos {\theta }_{\mathrm{2}}\sin {\alpha }_{\mathrm{2}}& -\cos {\theta }_{\mathrm{2}}\cos {\alpha }_{\mathrm{2}}& -\sin {\theta }_{\mathrm{3}}& \mathrm{1}\\ 0& 0& 1& 0\end{array}\right],\end{array}$$

(17)

where $r_i={(x_{1E}^2+y_{iN}^2+z_{iU}^2)}^{0.5}$ , ${\theta }_i$ and ${\alpha }_i$ are, respectively, the elevation and azimuth angle of the ith satellite. Since three satellites are in the same synchronization orbit, it is reasonable to assume that each satellite has the same pseudorange measurement error; namely, ${\sigma }_{\text{geo}}={\sigma }_1={\sigma }_2={\sigma }_3$ . Defining $w_h={\sigma }_h/{\sigma }_{\text{UERE}}$ , then W and the positioning accuracy ${\sigma }_p$ are, respectively, obtained as

$$\begin{array}{cc}& \mathbf{W}=\mathrm{diag}\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \mathrm{1}& w_h^2\end{array}\right],\\ {\sigma }_p=\text{WGDOP}\hspace{0.17em}{\sigma }_{\text{UERE}}=\sqrt{\mathrm{tr}{\left({\mathbf{H}}^T{\mathbf{W}}^{-\mathrm{1}}\mathbf{H}\right)}^{-\mathbf{1}}}{\sigma }_{\text{UERE}}.\end{array}$$

(18)

4.2. Performance Analysis of Emergency Positioning Algorithm with Altitude Observation

In terms of the integration of BDS with a baro-altimeter sensor aiding, different modeling errors are induced by applying different observation models to model the same altitude observation. Thus, the contribution to the positioning performance is also different. Therefore, in this section, we will focus on the contribution of the added altitude observation to positioning performance.

As the satellite pseudorange and aided observation have different measurement errors, it should be considered to analyze the contribution through weighting of the pseudorange and altitude measurements. In general, the position error can be expressed as

$$\begin{array}{l}{\mathbf{R}}_x={\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}}{\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{R}}_{\rho }{\mathbf{W}}^{-\mathrm{1}}\\ \hspace{1em}\times {\mathbf{H}}_N{\left({\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}}\right)}^T\\ ={\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}}{\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}\left[\begin{array}{ll}{\mathbf{I}}_{N-\mathrm{1}}{\sigma }_{\text{geo}}^{\mathrm{2}}& \mathbf{0}\\ \mathbf{0}& {\sigma }_h^{\mathrm{2}}\end{array}\right]\\ \times {\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N{\left({\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}}\right)}^T\\ ={\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}}{\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}\left(\mathbf{W}{{\sigma }_{\text{UERE}}^2}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}\right){\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\\ \times {\left({\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}}\right)}^T\\ ={\sigma }_{\text{UERE}}^{\mathrm{2}}{\left({\mathbf{H}}_N^T{\mathbf{W}}^{-\mathrm{1}}{\mathbf{H}}_N\right)}^{-\mathrm{1}},\end{array}$$

(19)

where N is the number of the observations for positioning. By substituting (18) into (19), this yields

$$\begin{array}{cc}& {\mathbf{P}}_N={\mathbf{P}}_N^T={\left({\mathbf{H}}_N^T{\mathbf{W}}^{-1}{\mathbf{H}}_N\right)}^{-1},\\ {\mathbf{P}}_{N-1}={\mathbf{P}}_{N-1}^T={\left({\mathbf{H}}_{N-1}^T{\mathbf{W}}_{N-1}^{-1}{\mathbf{H}}_{N-1}\right)}^{-1},\\ {\mathbf{h}}_N=\left[\begin{array}{cccc}{a}_{xN}& a_{yN}& a_{zN}& a_t\end{array}\right],\hspace{1em}\text{where}\hspace{0.17em}\hspace{0.17em}{a}_t=1\hspace{0.17em}\text{or}\hspace{0.17em}{a}_t=0.\end{array}$$

(20)

Since ${\mathbf{P}}_{N-1}$ is not singular, we simplify (19) by using the Sherman-Morrison formula, and we now have

$$\begin{array}{l}{\mathbf{R}}_x={\sigma }_{\text{UERE}}^2{\mathbf{P}}_N\\ ={\sigma }_{\text{UERE}}^2{\left({\mathbf{H}}_{N-1}^T{\mathbf{H}}_{N-1}\frac{{\sigma }_{\text{UERE}}^2}{{\sigma }_{\rho }^2}+{\mathbf{h}}_N^T{\mathbf{h}}_N\frac{{\sigma }_{\text{UERE}}^2}{{\sigma }_h^2}\right)}^{-1}\\ ={\sigma }_{\text{UERE}}^2{\mathbf{P}}_{N-1}-{\sigma }_{\text{UERE}}^2\frac{{\mathbf{P}}_{N-1}{\mathbf{h}}_N^T{\mathbf{h}}_N{\mathbf{P}}_{N-1}}{\left({\sigma }_h^2/{\sigma }_{\text{UERE}}^2\right)+{\mathbf{h}}_N{\mathbf{P}}_{N-1}{\mathbf{h}}_N^T}.\end{array}$$

(21)

According to (21), baro-altimeter observation's contribution $\mathrm{\Delta }C$ [22] can be defined as

$$\begin{array}{cc}& \mathrm{\Delta }C={\sigma }_{\text{UERE}}^2\text{trace}\left({\left[\frac{{\mathbf{P}}_{N-1}{\mathbf{h}}_N^T{\mathbf{h}}_N{\mathbf{P}}_{N-1}}{{\sigma }_h^2/{\sigma }_{\text{UERE}}^2+{\mathbf{h}}_N{\mathbf{P}}_{N-1}{\mathbf{h}}_N^T}\right]}_{\mathrm{3,3}}\right).\end{array}$$

(22)

As the elements in ${\mathit{h}}_N$ are not all 0, when $N>4$ , $\mathrm{\Delta }C>0$ [20], the larger the value, the more contribution to the positioning accuracy it has.

5. Experiment and Analysis

Through comparison and analysis, an improved altitude-aided station-centric-intersection emergency positioning algorithm is modeled. The below sections will further analyze and verify the conclusion above based on both simulation data and observation data. For convenience, the analysis assumes that the pseudorange and baro-altimeter altitude measurement errors are identically Gaussian distributed with zero-mean and variance ${\sigma }_{\text{geo}}^2$ , ${\sigma }_h^2$ [i.e., $N(0,{\sigma }_{\text{geo}}^2)$ ], respectively. Furthermore different measurements are statistically independent of each other. Meanwhile, to simplify the description below, three integrations are defined based on the five GEO satellites: 3GEOS-BAR, 4GEOS-BAR, and 5GEOS-BAR. Here, the 3GEOS-BAR is composed of G1, G3, G5, and baro-altimeter altitude sensor; 4GEOS-BAR is composed of G1, G3, G4, G5, and baro-altimeter altitude sensor; and 5GEOS-BAR is composed of five GEO satellites and baro-altimeter altitude sensor. The coordinate system is CGCS2000.

The BDS GEO observation data and baro-altimeter measurement data based on MS5607-02BA senor are collected by the ground station ST-A in China during 00:00 to 24:00 on August 7th, 2012, with a sampling rate of 1 Hz. Thus, there are 86400 epoch data to support analysis.

5.1. Raw Pseudorange Measurement Variance and Weighting Factor Selection

UERE is a measure of error in the range measurement to each satellite as seen by the receiver, which mainly contains control segment error, thermal noise, atmospheric error, and orbit mismatch error. As the difference of the satellite orbit determination error and extrapolation error of satellites in different orbit results in the unequal broadcast ephemeris accuracy, the UERE of GEO and IGSO will be unequal. To roughly estimate the satellite UERE, the observation data collected by ST-A is used for UERE statistics, the statistical formula can be expressed as [23]

$$\begin{array}{cc}& \text{UERE}={\rho }_{\text{cor}}-\left|{\mathbf{X}}_s-{\mathbf{X}}_u\right|,\end{array}$$

(23)

where ${\rho }_{\text{cor}}$ is the pseudorange corrected with satellite clock error, satellite relative group delay, ionospheric delay, troposphere delay, receiver clock error, and frequency deviation; ${\mathbf{X}}_s$ and ${\mathbf{X}}_u$ are the satellite and the receiver position vectors.

The statistics result based on regression cycle observation data are shown in Figure 6. It can be seen that the RMS of rough UERE (1σ) estimation for GEO and IGSO satellite are 6.5 m and 4.5 m, respectively. It should be noted that the UERE of IGSO satellite has ten-meter-level at 17:00, which is caused by broadcast ephemeris updating or satellite's entry into the service coverage area.

Figure 6

(a) The UERE based on raw observation data: GEO satellites. (b) The UERE based on raw observation data: IGSO satellites.

The principle of the baro-altimeter altitude measuring is to sense the atmospheric pressure so as to indirectly obtain the altitude. As the pressure and temperature vary with time, the altitude measurement error is deeply dependent on height, atmospheric pressure, and temperature, simultaneously. Statistic data shows that under the assumption that there is no temperature compensation and that the geodetic height is smaller than 1828.8 m, the relationship between the altitude h and the standard deviation ${\sigma }_h$ of the measured baro-altimeter error ${\epsilon }_h$ is as follows [21]:

$$\begin{array}{cc}& {\sigma }_h\left(h\right)=0.107\left(h-h_0\right)+{{\sigma }_h}_0\hspace{1em}\left(h<1828.8\hspace{0.17em}\text{m}\right),\\ {\sigma }_h\left(49\right)=457.2.\end{array}$$

(24)

Generally, the user altitude is below 100 m, thus, ${\sigma }_h$ will be smaller than 10.8 m, which can be obtained by substituting $h=100$  m into (24).

5.2. Positioning Algorithms Performance Analysis and Comparison Based on Different Altitude-Observation Models

In this section, in order to compare the performance of positioning algorithms using different altitude observation models, we will discuss the convergence characteristic and user positioning accuracy. For convenience, we will denote the positioning algorithms using the four altitude observation models above: approximate ellipsoid model, earth-centered fictitious satellite model, normal fictitious satellite model, and improved altitude observation model, as TM-1, TM-2, TM-3, and IM, respectively.

5.2.1. Analysis and Comparison of Convergence Characteristic

When Newton's iteration method is used to solve the equation set combined by altitude observation equation and pseudorange measurement equations, the initial position is required to be in close proximity to the user's true position, but this requirement is difficult to fulfill in practice. Therefore, reducing the requirement on the initial position is much more important. Under the condition of not changing the method, the requirement on the initial position can be reduced to some extent by using proper model. Thus, this section analyzes the effects of the initial position accuracy and the visible GEO satellite number on the convergence characteristic of the algorithm under the receiver cold's start condition.

Considering that the GEO satellite position deviates from the nominal equator orbit with some obliquity due to the various orbital perturbations during motion, the assumption that the simulation parameters, for example the number of visible satellite, the satellite position, and the user position are chosen based on that the first observation epoch data is more reasonable. And as noted in Section 5.1, the range errors, ${\sigma }_{\text{geo}}=6.5$ , ${\sigma }_h=10.8$ , would be used. The detail parameters are shown in Table 1.

Table 1

GEO satellite positions provided by broadcast ephemeris at the first epoch.

	Longitude (°)	Latitude (°)	H (m)
G1	139.99	1.59	35777978.92
G3	83.95	0.68	35788018.19
G4	159.99	0.26	35752467.17
G5	58.75	1.49	35789936.34

Two factors will be considered next. One is with a priori knowledge. The other is without altitude knowledge.

(1)  The Convergence Characteristic of Each Algorithm without Any A Priori Information. If no a priori knowledge is available, the initial position deviation $\mathrm{\Delta }\mathbf{x}$ needs to be searched in the 3-D space. Nonetheless, the region is limited in northern hemisphere to avoid converging to the mirror solution. Here, we use (25) to simulate the different deviation $\mathrm{\Delta }\mathbf{x}$ by changing k values. For simplification, Figure 7 illustrates each algorithm's convergence characteristic at $k=\mathrm{0.001}$ and $k=\mathrm{0.9}$ .

Figure 7

(a) Convergence performance comparison of different algorithms versus different initial positions: 3GEOS-BAR integration when $k=0.9$ . (b) Convergence performance comparison of different algorithms versus different initial positions: 3GEOS-BAR integration when $k=0.001$ . (c) Convergence performance comparison of different algorithms versus different initial positions: 4GEOS-BAR integration when $k=0.9$ . (d) Convergence performance comparison of different algorithms versus different initial positions: 4GEOS-BAR integration when $k=0.001$ .

Consider

$$\begin{array}{cc}& {\mathbf{x}}_0=k{\mathbf{x}}_u,\\ \parallel \mathrm{\Delta }\mathbf{x}\parallel =\parallel {\mathbf{x}}_c-\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{x}}_u\parallel ,\end{array}$$

(25)

where ${\mathbf{x}}_0$ , ${\mathbf{x}}_u$ , and ${\mathbf{x}}_c$ represent the initial position, the user's true position and the convergence result of the receiver position, respectively. $\mathrm{\Delta }\mathbf{x}$ represents the position deviation and k represents the deviation degree.

The convergence curve in Figure 7 shows that the convergence characteristic of different integrations is the same, independent to the integration. Generally, after 2 to 4 iteration, the algorithm has converged. When the initial position accuracy is high, all the algorithms can rapidly converge to a reasonable solution with the same accuracy. More specifically, when $k=0.001$ , although the deviation between the initial position and user's true position increases, all the algorithms successfully converge. However, only the convergence values of IM and TM-2 are the same as those of the situation $k=0.9$ , while TM-2 and TM-3 converge to unreasonable value. In terms of the convergence speed, there is no clear difference between these algorithms, all within 2–4 iteration. Table 2 shows the percentage of converging to reasonable solution when deviation degree k ranges from 0.001 to 100, with the variation step 0.001.

Table 2

Comparison of reasonable convergence area distribution without and with a priori knowledge, for different integration.

Algorithm	Sample size: $\mathit{3.6}e\mathit{5}$ sample points
	3GEOS-BAR		4GEOS-BAR
	Without a priori knowledge	With a priori known altitude	Without a priori knowledge	With a priori known altitude
IM	10.362%	49.154%	10.227%	49.154%
TM-1	6.169%	52.461%	6.457%	52.461%
TM-2	9.977%	48.122%	10.060%	48.122%
TM-3	10.358%	48.122%	10.222%	48.122%

When without any a priori knowledge, the initial position deviation $\mathrm{\Delta }\mathit{x}$ needs to be searched in the 3D space. Here, k is used to simulate the deviation distance in 3D direction. So the larger the k, the stronger the ability to resist the deviation. As shown in Figure 8(a) and Table 2, with the 3GEOS-BAR, the reasonable convergence range of IM, TM-2, and TM3 is more than 1.68 times of that of TM-1. More specifically, the reasonable range of TM-1 is only $0.15\leqq k\leqq 6.303$ , corresponding to $0\leqq \parallel \mathrm{\Delta }\mathbf{x}\parallel <3.37e7$  m. while the range of other algorithms’ convergence region is $0.08<k<10.5$ , corresponding to $0\leqq \parallel \mathrm{\Delta }\mathbf{x}\parallel <6.05e7$  m. In terms of integrations, by adding one GEO satellite into the integration, slightly negative effects would be caused on IM and TM-3.

Figure 8

(a) Reasonable convergence area distribution with the 3GEO-BAR integration: without any a priori knowledge. (b) Reasonable convergence area distribution with the 3GEO-BAR integration: with a priori known altitude using IM. (c) Reasonable convergence area distribution with the 3GEO-BAR integration: with a priori known altitude using TM-2.

(2)  The Convergence Characteristic of Each Algorithm with Altitude A Priori Information. When with a priori knowledge, obtained from an external sensor, the search range would be limited in the horizontal plane by using ENU coordinates to analyze the convergence characteristic accordingly. Thus, the initial position deviation $\mathrm{\Delta }\mathbf{x}$ can be defined as (26)

$$\begin{array}{cc}& x_{\mathrm{0}}=R\cos \left(\theta \right),\hspace{1em}\mathrm{0}\le \theta \le \mathrm{2}\pi ,\\ y_{\mathrm{0}}=R\sin \left(\theta \right),\\ z_{\mathrm{0}}=\mathrm{0},\end{array}$$

(26)

where $x_0$ , $y_0$ , and $z_0$ are the 3D components of the initial position vector ${\mathbf{x}}_{\mathrm{0}}$ . θ is the angle between the 3D vector direction and the east direction in ENU coordinates, and $R=\parallel \mathrm{\Delta }\mathbf{x}\parallel$ . In order to compare with the result in Section 5.2.1(1), the upper limit of R is adjusted to be larger than the maximum convergence range ${\parallel \mathrm{\Delta }\mathbf{x}\parallel }_{\max }=6.05e7$  m in Section 5.2.1(1). Thus, the following parameters are used: $0<R\leqq 1e8$  m, $\mathrm{\Delta }R=1e4$  m, $0<\theta \leqq {360}^{°}$ , and $\mathrm{\Delta }\theta =1^{°}$ .

After comparing the ability to resist the initial position deviation between the conventional algorithms and the improved algorithm under the same condition, it can be seen that the larger the reasonable convergence region area, the stronger the ability to resist the deviation. The detailed result is shown in Table 2. More specifically, with the integration of 3GEOS-BAR, in comparison to the three conventional algorithms, the convergence region of IM is, respectively, increased by −3.307%, 1.03%, 1.03%, and which is weaker than that of TM-1. However, when adding a GEO satellite into the 3GEOS-BAR integration, it will result in the convergence area of TM-1 decreasing by 44.26%, while the convergence region of the other three algorithms respectively increasing by 5.82%, 6.68%, 6.73%. Therefore, in comparison to TM-2 and TM-3, IM has some advantages, but this advantage weakens with the satellite number increasing, which complies with the conclusion in Section 5.2.1(1). Nonetheless, the antideviation ability of TM-1 rapidly decreases with the satellite number increasing. In other words, the robustness of TM-1 is very poor, which may limit its application scope.

Then we would like to analyze the reasonable convergence region distribution within a plane, which is shown in Figures 8(b) and 8(c). For simplification, only the analysis of Figure 8(b) would be performed next. And the analysis of Figure 8(c) is similar.

Referring to Figure 8(b). The reasonable convergence region is not coordinate system symmetrical, while it is impacted by R and θ simultaneously and has obvious directivity. Under the same condition, the North direction is more sensitive to the R and θ variation, which results in the area ratio of the two reaching to 1 : 6.44; and with R gradually approaching to its limit, the ratio will continuously increase. Further simulation shows that only when R reaches $9.06e12$  m, it will converge to unreasonable value. In comparison to the situation without a priori knowledge in Section 5.2.1(1), the ability to resist the initial position deviation of each algorithm is evidently increased.

5.2.2. Analysis and Comparison of Positioning Accuracy

To analyze the effects of the modeling error of different altitude observation models on the positioning accuracy, the two related factors of the user receiver altitude and integration would be considered in the simulations. Here, the orbit data is generated by BDS GEO standard constellation, pseudorange is the value of geometrical distance from the receiver to each satellite adding with the pseudorange measurement error, and altitude observation value is the value of true altitude adding with the baro-altimeter altitude measurement error; Meanwhile, to avoid the initial position deviation $\mathrm{\Delta }\mathbf{x}$ being too large to result in the positioning solution converging to an unreasonable value, the $\parallel \mathrm{\Delta }\mathbf{x}\parallel \leqq 100$  m is chosen.

(1)  Each Algorithm's Positioning Accuracy at Different Altitude. As an example of the positioning accuracy at different altitude, consider the 3GEOS-BAR. In the simulation, the altitude of the user (46°N, 105°E) ranges from zero to 1.2 times of the GEO satellite altitude, and the variation step is $1e4$  m. At each altitude value, $1e5$ sample points are generated for statistical analysis. Here, the mean and standard deviation of positioning error are obtained. To demonstrate the accuracy difference of different algorithms within the entire altitude range more accurately, the simulated positioning result data is divided into two groups according to the altitude at the rate of $1e6$  m: the first group with 25 sample points, corresponding to $0\leqq h\leqq 2.5e7$ ; and the second group with 18 sample points, corresponding to $2.5e7<h\leqq 4.2e7$ . The simulation result is shown in Figure 9.

Figure 9

(a) Comparison of the mean and standard deviation of the positioning error, for different altitude observation models: mean value when $0\leqq h\leqq 2.5e7$ . (b) Comparison of the mean and standard deviation of the positioning error, for different altitude observation models: mean value when $2.5e7<h\leqq 4.2e7$ . (c) Comparison of the mean and standard deviation of the positioning error, for different altitude observation models: standard deviation when $0\leqq h\leqq 2.5e7$ . (d) Comparison of the mean and standard deviation of the positioning error, for different altitude observation models: standard deviation when $2.5e7<h\leqq 4.2e7$ .

Figures 9(a) and 9(b) show that IM can always compute the normal position fix, and the mean error is close to zero, which is not affected by the altitude variation. While the mean error of TM-1 is seriously affected by the user altitude. When the altitude reaches $3.58e7$  m, the direction cosine matrix H of the algorithm becomes singular to cause failures of computing the position fix (marked with light green boxes in the bottom of Figure 9(b)). And the performance of TM-2 and TM-3 is the same, but different from IM. Specifically, when the altitude reaches $1.9e7$  m, deviation appears simultaneously at the East-West direction and North-South direction, particularly, and the deviation at North-South direction increases more quickly, leading to the wrong positioning result.

Further, Figures 9(c) and 9(d) depict the standard deviation variation of the positioning error. The error in the altitude direction of IM is the same as the chosen value ${\sigma }_h$ , and it is independent of the altitude variation, while the errors in East-West direction and North-South direction decrease and then increase with the altitude varying. Note that the error at North-South direction is larger than that at East-West direction. Nonetheless, the standard deviation variation of TM-1 is not the same as that of mean error, only the vertical error is caused. The difference of TM-2, TM-3, and IM is similar to that of the mean positioning error. The reason is that the geometry between the satellites and the user varies with the altitude continuously increasing, which results in the positioning error also varying accordingly. And the reason for the standard deviation in the North-South direction being much larger than that in East-West direction is mainly that the GEO satellites in BDS are all located at the equator with the same latitude, which makes North dilution of precision (NDOP) much larger than East dilution of precision (EDOP).

Further, the performance of different algorithm is shown in Table 3. In addition, the availability within the entire altitude range is also provided, which is obtained under the assumption that the horizontal positioning error threshold is 65 m based on the analysis and statistics of HDOP and ${\sigma }_{\text{UERE}}$ in other sections.

Table 3

Performance comparison of different positioning algorithms.

Algorithm	RMS error (m)				Availability
	East	North	UP	3D
IM	8.479	43.335	10.858	45.472	93.3%
TM-1	8.440	45.024	211.883	216.787	77.1%
TM-2	8.478	48.863	10.859	50.768	61.4%
TM-3	8.478	48.863	10.859	50.768	61.4%

It can be seen from data in Table 3 that, in comparison to conventional TM-1, TM-2, and TM-3, the improvements in 3D positioning accuracy is 79.0%, 10.4%, and 10.4%, respectively, and improvements in the availability is 16.2%, 15.7%, and 15.7%, respectively. Thus, the performance advantage is obvious. It should be noted that the performance of TM-2 and TM-3 is very close in each, with the maximum difference of the error below $1e-6$  m.

(2)  Positioning Error Contribution of Each Algorithm with Different Integration. In the simulations, the visible GEO satellite number for different integrations would range from 3 to 5, corresponding to 3GEOS-BAR, 4GEOS-BAR, and 5GEOS-BAR integrations, respectively, and user's position is chosen at (46°N, 105°E, and 100 m). For each integration, $1e5$ sample points of simulation data are generated, and each sample point includes the pseudorange and altitude data whose quantity is the same as that of visible satellites. The root-mean-squared (RMS) error statistics in each direction is shown in Figure 10.

Figure 10

(a) Comparison of the RMS error for different altitude observation models: East. (b) Comparison of the RMS error for different altitude observation models: North. (c) Comparison of the RMS error for different altitude observation models: Vertical. (d) Comparison of the RMS error for different altitude observation models: 3D.

Figure 10 shows that for 3GEOS-BAR integration, the errors of all algorithms are very close, and the difference is only at $1e5$  m level. When adding G4 satellite into the 3GEOS-BAR integration, the positioning accuracy of IM, TM-1, TM-2, and TM-3 is improved by different extents which are 82.20%, 56.09%, 82.19%, and 82.20%, respectively. Furthermore, if G2 is further added into the integration, the four values are, respectively, 48.85%, 31.302%, 48.85%, and 48.85%. More specifically, the errors of IM and TM-3 are always the smallest, and the error of TM-2 is slightly larger than that of two, which is only about $1e-2$ meter. However, the errors in all directions of TM-1 are the largest, and only the error in the East-West direction is basically close to that of the other two algorithms.

Furthermore, according to (22) in Section 4.2, the theoretical 3D position error contribution can be obtained. With the simulation result above, the comparison results are shown in Table 4, where the subscript i represents the integration with the i GEO satellite. $\mathrm{\Delta }{C}_{i,j}^T$ or $\mathrm{\Delta }{C}_{i,j}^{\text{rms}}$ ( $3⩽i⩽5$ , $3⩽j⩽5$ , $j<i$ ) represents the 3D position error contribution of integration i relative to the integration j, and the superscript rms and T represent the statistics value and theoretical value, respectively.

Table 4

Comparison of different altitude observation model contributions to position error obtained from simulation and theoretical results, for different algorithms, for different integration.

Algorithm	$\mathrm{\Delta }{C}_{\mathrm{4,3}}^{\text{rms}}$	$\mathrm{\Delta }{C}_{\mathrm{4,3}}^T$	$\mathrm{\Delta }{C}_{\mathrm{5,4}}^{\text{rms}}$	$\mathrm{\Delta }{C}_{\mathrm{5,4}}^T$
IM	236.958	221.835	25.063	31.326
TM-1	161.692	161.174	39.621	40.287
TM-2	236.940	227.713	25.073	31.534
TM-3	236.958	222.193	25.063	30.967

From Table 4, it can be seen that the position error contribution obtained from statistics complies with the theoretical value, thus, verifying the correctness of the theory analysis in this paper. It also states that under the same condition, the positioning error contribution of TM-1 is the least.

5.3. Analysis of Satellite Visibility with Baro-Altimeter Sensor Aiding

At present, the provided RNSS and RDSS service coverage range based on GEO satellites is for China and the surrounding area [19], specifically, between 5°N~55°N and 70°E~140°E [23]. For the 3D passive positioning, four satellite coverage is at least required. While, in principle, the baro-altimeter sensor aiding is equivalent to adding a fictitious satellite at the Earth center, and then the visibility will change. Therefore, to determine the variation of visibility, we ran a simulation with different observation locations in a 0.5° grid from 90°N to 90°S and 180E°–180°W and in 1 minute increments over a day, with 5° and 15° mask angle, respectively. Take 5° mask elevation angle as an example, Figure 11(a) and Table 5, show the 5GEOS-BAR visibility and the statistics result, respectively.

Table 5

Visibility of 5GEOS-BAR integration with 5° and 15° mask angles.

Integration	Mask angle (°)	Percentage of at least N visible satellite
		≥3	≥4	≥5
5GEOS-BAR	5	32.4%	18.3%	9.2%
	15	23.3%	11.4%	4.1%

Figure 11

(a) The average number of visible satellites with 5 GEO satellites. (b) Cumulative distribution of HDOP with and without altitude-aiding.

It can be seen that the percentage of three, four, and five satellite coverages is, respectively, 32.4%, 18.3%, and 9.2%. Thus, five satellite coverages for the China region can be realized. However, the latitude range for three satellite coverage is only 76°S~76°N, which leads to a low-coverage for high latitude area. Figure 11(b) illustrates HDOP cumulative distribution of HDOP with and without baro-altimeter sensor aiding. The global availability in the three cases is all very low. If the HDOP threshold 10 is chosen, with the baro-altimeter sensor aiding, the availability at 5° and 15° mask angle is 9.81% and 6.3%, respectively. If baro-altimeter sensor aiding is removed, the availability is reduced to 4.42%. Therefore, obvious availability improvement of the GEO satellite emergence service through baro-altimeter sensor aiding is realized.

5.4. Performance Comparison between Altitude-Aided Station-Centric-Intersection Emergency Positioning Algorithm and Conventional Algorithm

In comparison to the conventional positioning algorithm, the difference of the improved altitude-aided station-centric-intersection emergency positioning algorithm is through weighting of the baro-altimeter measurement in ENU coordinates to improve user's positioning accuracy. Here, to verify the performance of the improved algorithm, the method of changing the weighting factor $w_h$ is applied for different integrations. At each $w_h$ , $1e5$ sample points of simulation data are generated, the other simulation parameters are the same as those in Section 5.2.2. According to the navigation environment and baro-altimeter accuracy specification, the variation range of $w_h$ is from 0.1 to 20, and the step is 0.1. Then, comparison of the positioning RMS error between the improved algorithm and the conventional algorithms is performed, which is shown in Figure 12, where $E_W$ and $E_T$ represent the RMS error of the improved algorithm and the conventional algorithm, respectively. More specifically, it can be seen that when the number of visible GEO satellites is more than 3, the RMS error difference $E_W-E_T$ increases with $w_h$ . In other words, the improvement of positioning accuracy is remarkable when the baro-altimeter sensor measurement error greatly differs from the GEO pseudorange measurement error. However, when the number of visible GEO satellites is only 3, there is no obvious difference of the two algorithms.

Figure 12

(a) Comparison of positioning error difference between the improved algorithm and conventional algorithms of different integration with simulation data: East. (b) Comparison of positioning error difference between the improved algorithm and conventional algorithms of different integration with simulation data: North. (c) Comparison of positioning error difference between the improved algorithm and conventional algorithms of different integration with simulation data: Vertical. (d) Comparison of positioning error difference between the improved algorithm and conventional algorithms of different integration with simulation data: 3D.

Therefore, in various applications, the improved algorithm should be adopted when a higher reliability is required on the positioning accuracy, such as aircraft landing, which enhances accuracy through weighting of pseudorange and baro-altimeter measurements based on their relative noise levels and relative importance for each estimated quantity.

5.5. Experimental Results

In order to confirm the conclusions above, which are obtained based on simulation pseudorange data, experimental analysis is performed, based on 24 hours of BDS observation data collected by ST-A. Here, limited by the length of the paper, we will focus on the effects of different altitude observation models on positioning accuracy and performance comparison between the altitude-aided station-centric-intersection emergency positioning algorithm and the conventional algorithm.

(1)  Effects of Different Altitude Observation Model on Positioning Accuracy. To avoid the initial position deviation $\mathrm{\Delta }\mathbf{x}$ being too large to result in the positioning solution converging to an unreasonable value, the $k=0.9$ in (25) is chosen. Figure 13 shows the 3D RMS error of different algorithms with the 4GEOS-BAR integration, where the x-axis represents the Beidou navigation satellite system Time (BDT) seconds of week. Figure 13(a) shows the 3D RMS errors of IM are on the order of 5.909 to 74.727 m, and the points representing the error of TM-2 and TM-3 almost overlap with those of IM. In other words, the difference among IM, TM-2, and TM-3 is very small. However, for TM-1, many positioning failure occurs, as shown in Figure 13(b). In order to evaluate the performance of TM-2 completely, the error threshold $V_t=200$  m is applied to IM-2 result, and any 3-D RMS error that is at or above the $V_t$ is removed; hence, the statistics result based on the left data is obtained, as shown in Table 6.

Table 6

Comparison of the positioning errors and failure rate of positioning for different algorithms with the 4GEOS-BAR integration.

Algorithm	RMS error (m)				Failure rate of positioning
	East	North	UP	3-D
IM	6.343	14.863	10.826	19.451	0
TM-1	6.843	69.900	63.201	94.484	18.62%
TM-2	6.343	14.864	10.826	19.452	0
TM-3	6.343	14.864	10.826	19.452	0

Figure 13

(a) Comparison of the 3D RMS error for different algorithms with the 4GEOS-BAR integration: IM, TM-2, and TM-3. (b) Comparison of the 3D RMS error for different algorithms with the 4GEOS-BAR integration: TM-1.

Table 6 shows that when the position fix is based on 4GEOS-BAR integration, RMS errors in East, North, vertical, and 3D direction of IM, TM-2, and TM-3 are almost the same. However, only the East RMS error of IM-1 is the same as those of the three algorithms, and the errors in other directions remarkably increase to about 5 times. This result complies with the one in Section 5.2.2; hence, the correctness of the simulation conclusion is verified.

(2)  Performance Comparison between the Altitude-Aided Station-Centric-Intersection Emergency Positioning Algorithm and the Conventional Algorithm. According to the simulation method in Section 5.4, Figures 14(a) and 14(b) demonstrate the positioning error of the altitude-aided station-centric-intersection emergency positioning algorithm (denoted as improved algorithm in Figure 14) and the conventional algorithm with the two integrations of 3GEOS-BAR and 4GEOS-BAR, based on BDS observation data, respectively. For simplification, only the 3D direction performance comparison is illustrated.

Figure 14

(a) Comparison of the 3D position errors between the improved algorithm and the conventional algorithms with two integrations, based on with observation data: 3GEOS-BAR integration. (b) Comparison of the 3D position errors between the improved algorithm and the conventional algorithms with two integrations, based on with observation data: 4GEOS-BAR integration.

As shown in Figure 14, with the integration of 3GEOS-BAR, the 3D position error difference of the two algorithm positioning errors is not obvious. However, when a GEO satellite is added into the integration, in comparison to the conventional algorithm, the new altitude-aided station-centric-intersection emergency positioning algorithm can greatly improve the positioning accuracy of the solution, and the greater the difference between the satellite pseudorange and baro-altimeter altitude measurement, the more significant the effect. Taking the $w_h=15$ with the 4GEOS-BAR integration as an example, the positioning error of the conventional algorithm is 164.897 m while the one of the improved algorithm is only 93.205 m, and 43.477% improvement is achieved. Hence, the reliability of evaluating the positioning accuracy is also greatly improved. This result complies with the simulation result in Section 5.2.2, so the correctness of the simulation result is verified.

6. Conclusion

BDS uses five GEO satellites to provide RNSS and RDSS service simultaneously, so as to take the advantages of fixed coverage area and high navigation message speed of GEO satellites. In practice, the poor user/GEOs relative geometry results in severe user position error, even unavailability. If an external sensor can provide altitude information, 2D BDS positioning can be a remedy for some application, providing an improvement in system performance. However, due to the modeling error of the conventional baro-altimeter observation models, the positioning accuracy and availability cannot completely meet the requirement in some cases, such as the user's altitude exceeding over $1e7$ meter and more than 3 GEO satellites in the integration. This paper analyzes the principle of GEO satellites and baro-altimeter sensor integration and the effects of the standard deviation of the baro-altimeter measurement error on positioning errors, and proposes an altitude-aided station-centric-intersection emergency positioning algorithm through weighting of pseudorange and altitude measurements. The algorithm is accurate, without modeling error, and is not limited by user's altitude and the integrations. By simulation and experiment, it is shown that the proposed algorithm can effectively enhance the positioning accuracy, availability, and robustness.

Considering the launch time of GEO-5 and the observation data collected time in this paper, there is no GEO-5's observation data for analysis in this paper. So future work is being carried out to validate the performance of the proposed algorithm with the 5GEOS-BAR integration, based on BDS observation data. Also the advantages that IGSO satellites carry RDSS payload are being analyzed.