*2.3. DGNSS/DBA Combined Positioning Algorithm*

In this study, we propose a DGNSS/DBA combined positioning algorithm, in which the user's altitude obtained by the DBA system in Section 2.2 is used to construct the Earth ellipsoid constraint equation and solved rigorously as an independent observation jointly with the DGNSS observation equation, which is equivalent to adding a virtual satellite located at the center of the Earth [26]. Since the geodetic height is an independent variable in GNSS coordinates, an approximate ellipsoid with the altitude *h* from the reference ellipsoid (WGS-84) can be constructed using the user geodetic height, as shown in Figure 1.

**Figure 1.** Approximate reference ellipsoidal meridian profile where the ground user's geodetic height is located.

At this time, when the ground user's geodetic height is not very large, the observation equation after DGNSS/DBA combination can be expressed as:

$$\begin{cases} \nabla \Delta P\_{i,j}^{s\_1 s\_k} = \nabla \Delta \rho\_{i,j}^{s\_1 s\_k} + \nabla \Delta \varepsilon\_{i,j}^{s\_1 s\_k} \\ \frac{\lambda^2 + \lambda^2}{\left(a + h\right)^2} + \frac{Z^2}{\left(b + h\right)^2} = 1 \end{cases} \tag{4}$$

the symbols in the DGNSS observation equation in the first line of Equation (4) are the same as Equation (2). *P*(*X*,*Y*, *Z*) is the 3D coordinate of the ground user; *a* and *b* are the long and short semiaxes of the WGS-84 Earth reference ellipsoid, respectively. Since *h* is much smaller than the long and short semiaxes of the Earth reference ellipsoid, an approximate reference ellipsoid with a long semiaxis *a* + *h* and short semiaxis *b* + *h* is used instead without causing much bias [37]. To solve the Earth ellipsoid constraint equation in the second expression of Equation (4) by differential processing, the ellipsoid constraint equation is expanded in the user's approximate position (*X*0,*Y*0, *Z*0) according to the Taylor series, and only the first-order term is retained, where the partial derivative of *X* is obtained as:

$$\frac{2X\_0}{\left(a+h\right)^2}dX - \frac{2X\_0^2}{\left(a+h\right)^3}dh - \frac{2Y\_0^2}{\left(a+h\right)^3}dh - \frac{2Z\_0^2}{\left(b+h\right)^3}dh = 0\tag{5}$$

after simplification, we get:

$$\frac{\partial h}{\partial X\_0} = \frac{X\_0(a+h)(b+h)^3}{\left(X\_0^2 + Y\_0^2\right)\left(b+h\right)^3 + Z\_0^2\left(a+h\right)^3} \tag{6}$$

similarly, taking partial derivatives of *Y* and *Z* yields:

$$\frac{\partial h}{\partial Y\_0} = \frac{Y\_0(a+h)(b+h)^3}{\left(X\_0^2 + Y\_0^2\right)(b+h)^3 + Z\_0^2(a+h)^3} \tag{7}$$

$$\frac{\partial h}{\partial Z\_0} = \frac{Z\_0(a+h)^3(b+h)}{\left(X\_0^2 + Y\_0^2\right)\left(b+h\right)^3 + Z\_0^2\left(a+h\right)^3} \tag{8}$$

let *α* = *∂h*/*∂X*0, *β* = *∂h*/*∂Y*0, *γ* = *∂h*/*∂Z*0, that is, the Earth ellipsoidal constraint equation in Equation (4) is linearized at the approximate coordinates (*X*0,*Y*0, *Z*0) to give:

$$V\_{DBA} = \alpha dX + \beta dY + \gamma dZ - d\hbar\tag{9}$$

where *dh* <sup>=</sup> *<sup>h</sup>* <sup>−</sup> <sup>ˆ</sup> *h* is the altitude residual, *h* is the altitude obtained by the DBA system, and ˆ *h* is the geodetic height obtained by the users' at the approximate position. *dX*, *dY*, *dZ* are positional corrections of the receiver antenna center. The detailed conversion process can be found in the original literature [38].

Similarly, the DGNSS observation equation in the first expression of Equation (4) is expanded by the Taylor series at the approximate position (*X*0,*Y*0, *Z*0), omitting higher-order terms above the first order, and combined with Equation (9) to obtain the DGNSS/DBA combined positioning error equation:

$$V = H\mathfrak{X} - l,\\
P\tag{10}$$

in Equation (10), the parameter estimated as *x*ˆ = *dX dY dZ <sup>T</sup>* contains three approximate position correction values; *<sup>V</sup>* <sup>=</sup> *<sup>V</sup>*<sup>1</sup> ··· *<sup>V</sup><sup>i</sup> VDBA <sup>T</sup>* is the residual vector, ⎡ <sup>1</sup> *m*<sup>1</sup> *n*<sup>1</sup> ⎤

*H* = ⎢ ⎢ ⎣ *l* ··· ··· ··· *l <sup>i</sup> m<sup>i</sup> n<sup>i</sup> αβγ* ⎥ ⎥ ⎦ is the coefficient matrix, *l <sup>i</sup>* = (*Xs*−*X*0) *ρ* (0)*s j* <sup>−</sup> (*Xk*−*X*0) *ρ* (0)*k j* , *<sup>m</sup><sup>i</sup>* <sup>=</sup> (*Ys*−*Y*0) *ρ* (0)*s j* − (*Yk*−*Y*0) *ρ* (0)*k j* , *<sup>n</sup><sup>i</sup>* <sup>=</sup> (*Zs*−*Z*0) *ρ* (0)*s j* <sup>−</sup> (*Zk*−*Z*0) *ρ* (0)*k j* are the pseudorange double-difference directional cosine, respectively. *<sup>l</sup>* <sup>=</sup> *<sup>L</sup>*<sup>1</sup> ··· *<sup>L</sup><sup>i</sup> dh <sup>T</sup>* are the observation value vectors. *<sup>P</sup>* <sup>=</sup> *PDGNSS* <sup>0</sup> <sup>0</sup> *PDBA* is the DGNSS/DBA combined positioning weight matrix. *PDGNSS* = *Q*−<sup>1</sup> *DD* is the a priori weight matrix of GNSS pseudorange double-difference observation

equation. *QDD* is the GNSS pseudorange double-difference observation values covariance matrix, the stochastic model of GNSS nondifferential observations adopts the sine trigonometric function elevation angle fixed-weight model [39]. According to the error propagation law, the GNSS relative positioning variance-covariance matrix can be expressed as *QDD* [34]. The a priori weight *PDBA* of the DBA can be determined based on the results of the empirical evaluation in Section 3.2. We can solve Equation (10) using the single-epoch weighted least squares method, as Equation (11):

$$\begin{cases} \text{ } \pounds = \left(H^T P H\right)^{-1} H^T P l\\ \quad Q\_{\pounds} = \left(H^T P H\right)^{-1} \end{cases} \tag{11}$$

where *Qx*<sup>ˆ</sup>*x*<sup>ˆ</sup> is the posterior covariance matrix of the parameter *x*ˆ. It can be found that the DBA altitude constraint is equivalent to adding a code pseudorange observation, and the condition number of the error equation coefficient matrix *H* becomes significantly better, the position dilution of precision (PDOP) value can be effectively reduced. However, the accuracy of the positioning solution is influenced by the accuracy of DBA altitude, i.e., if the accuracy of DBA altitude is better than DGNSS altitude, the improvement effect is obvious, otherwise, the positioning accuracy cannot be improved.

#### **3. Experiment Results**

#### *3.1. The Introduction of Experiment Data*

In the experiment, the reference station included a high-precision geodetic GNSS receiver Trimble NET R9, a Trimble Choking 59,800 antenna, and a low-cost BMP280 barometer. The mobile station consisted of a single-frequency low-cost u-blox receiver NEO-M8T, a patch antenna, and two BMP280 barometers. The mobile station was also equipped with a geodetic receiver Trimble NET R9 for data acquisition, and the postprocessed kinematic (PPK) mode of commercial software Inertial Explorer 8.70 was used to process data to obtain high-precision 3D coordinate sequences as reference values. The nominal resolution of the low-cost BMP280 barometer was 0.01 mbar (0.1 m) and the data sampling rate was set at 1 Hz. As shown in Figure 2, to prevent the effect of crosswind on the barometric pressure, the barometer was placed inside a transparent plastic box with several small holes at the top of the box. All experimental data were recorded and postprocessed by a laptop computer. The height deviation of the BMP280 barometer from the GNSS antenna phase center was compensated by data preprocessing, and the DGNSS/DBA combined positioning analysis was performed by the self-written programs.

The whole experiment was divided into three parts. The first experiment was used to evaluate the DBA altitude accuracy at different baseline lengths and to provide a priori information for subsequently combined positioning with DGNSS. The second experiment evaluated the positioning performance of low-cost single-frequency DGNSS and DGNSS/DBA combined positioning through static experiments with 65 m and 6.0 km baseline lengths. The third experiment evaluated the kinematic vehicle positioning performance of low-cost single-frequency DGNSS and DGNSS/DBA combined positioning in open and complex urban environments, respectively.
