A Novel All-Weather Method to Determine Deflection of the Vertical by Combining 3D Laser Tracking Free-Fall and Multi-GNSS Baselines

Abstract

The bright stars in the clear night sky with weak background lights should be observed in the traditional deflection of the vertical (DOV) measurement so that the DOV cannot be observed under all-weather conditions, which limits its wide applications. An all-weather DOV measurement method combining three-dimensional (3D) laser tracking free-fall and multi-GNSS baselines is proposed in this paper. In a vacuum environment, the 3D laser tracking technique is used to continuously track and observe the motion of free-fall with high frequency and precision for obtaining 3D coordinate series. The plumb line vector equation is established to solve the gravity direction vector in the coordinate system of the laser tracker at the measuring point using least squares fitting coordinate series. Multi-GNSS observations are solved for obtaining the precise geodetic cartesian coordinates of the measuring point and GNSS baseline information. A direction transformation method based on the baseline information proposed in this paper is used to convert the gravitational direction vector in the laser tracker coordinate system into the geodetic cartesian coordinate system. The geodetic cartesian coordinates of the measuring point are used to calculate the ellipsoid normal vector, and the angle between this and the gravity direction vector in the geodetic cartesian coordinate system is estimated to obtain the astrogeodetic DOV. The DOV is projected to the meridian and prime vertical planes to obtain the meridian and prime vertical components of the DOV, respectively. The astronomical latitude and longitude of the measuring point are calculated from these two components. The simulation experiments were carried out using the proposed method, and it was found that the theoretical precision of the DOV measured by the method could reach 0.2″, which could realise all-weather observation.

1. Introduction

The deflection of the vertical (DOV) is the angle between the plumb line and the normal to the reference ellipsoid, indicating the inclination of the geoid with respect to the reference ellipsoid, and characterising the spatial and temporal inhomogeneity of the mass distribution within the earth [1,2]. The DOV is indeed a basic observation of geodesy [3,4,5], which determines the gravity direction and contains rich high-frequency information about the gravity field [6]. It can be used to infer the size and shape of the mean earth ellipsoid [7,8], refine the geoid model, convert astrogeodetic data [9,10], and correct precision engineering measurements [11]. The time-varying study of DOV verifies that there is a very close relationship between its change and earthquakes, which provides a new means of monitoring large earthquakes [12,13]. In conclusion, the DOV can provide rich detailed information about the earth’s gravity field and the geoid, which is of great significance in the fields of geodetic applications, geophysical inversion, resource exploration, and seismic and volcano monitoring [2,14,15,16].

The common methods for measuring DOV include astronomical geodesy, gravimetry, astronomical gravimetry, and GNSS levelling [1,17]. Astronomical geodesy generally uses high-precision astronomical theodolites or digital zenith camera systems [18,19,20] to observe low-magnitude stars for achieving high-precision DOV [21,22]. However, such instruments are cumbersome and inefficient in measurement, and can only observe lower-magnitude stars on clear nights, not all-weather. The gravimetry essentially uses methods such as the Stokes formula or the Vening Meinesz formula to integrate gravity anomalies on geoids to calculate the DOV [23,24]. However, the method assumes no disturbances outside the geoid and global gravity anomalies are known, neither of which can be achieved, so it is rarely used independently. The astronomical gravimetry method is a combination of astronomical geodesy and gravimetry to determine DOV. Therefore, this method, like astronomical geodesy, is seriously constrained by the measurement environment. The GNSS-levelling method uses height anomaly to calculate DOV [1,25], requiring the survey area to be flat. The measurement method is relatively low in accuracy and efficiency, which has some limitations. In general, these traditional methods to determine DOV have problems such as difficulty in measurement, being time consuming, susceptibility to environmental conditions, and an inability to observe under all-weather conditions. It seriously limits the application of DOV. Therefore, all-weather and high-precision DOV measurement is still a pressing challenge in the world [2,15,22].

The three-dimensional (3D) laser tracking technique is a polar coordinate-measuring method with high measuring speed and precision. With characteristics such as a data sampling rate of up to 1000~3000 Hz, the laser tracker measurement system with sub-micrometre level precision can track the target automatically and realise accurate and dynamic three-dimensional tracking measurement [26,27]. Based on this, an all-weather DOV measurement method combining 3D laser tracking free-fall and multi-GNSS baselines is proposed in this paper.

The structure of this paper is as follows: Section 2 introduces the principles for the all-weather method to determine DOV by combining 3D laser tracking free-fall and multi-GNSS baselines. Section 3 conducts simulation experiments to analyse the theoretical precision of the proposed method, and separately discusses the influence of laser tracking and GNSS errors on DOV. Section 4 contains the conclusion.

2. Principle and Methodology

The schematic diagram and technical route of DOV measurement using the proposed method are shown in Figure 1 and Figure 2. The method involves two kinds of coordinate systems: one is the spatial cartesian coordinate system established with the centre of the laser tracker as the origin, called laser tracker coordinate system (LTCS), and the other is the geodetic cartesian coordinate system (GCCS). Firstly, the central coordinates of the GNSS antennas are measured using the laser tracker, and they are differentiated to obtain the baseline information in LTCS. The laser tracker is used to dynamically track and measure the 3D coordinates of the free-fall target. Then, the gravity direction vector in LTCS is calculated using least squares (LS) fitting 3D coordinate series. With the sub-micrometre precision distance constraint of baseline, the GNSS data are processed to obtain high-precision geodetic coordinates (or geodetic cartesian coordinates) of the measuring point P and the GNSS antennas’ baseline vectors in GCCS. According to the baseline information in LTCS, the direction transformation parameters between LTCS and GCCS are calculated to convert the gravity direction vector in LTCS into GCCS. According to the ellipsoid normal vector determined by the geodetic coordinates of the point P, the astrogeodetic DOV is estimated and projected onto the meridian and prime vertical planes to obtain the meridional and prime vertical components of the DOV, respectively. Finally, the astronomical coordinates are obtained from the two components and the geodesic coordinates of the point P.

2.1. Extracting Gravity Direction Vector

In a vacuum environment, the laser tracker is used to track the free-falling target with high precision and dynamics for obtaining its three-dimensional coordinate series in LTCS. According to the coordinate series, the LS and total least squares (TLS) fitting methods [28,29,30] is used to calculate the gravity direction vector in LTCS, respectively.

2.1.1. Spatial Linear Fitting Using LS

Assume that the coordinates of the free-falling target in LTCS are $L_i\left(x_i,y_i,z_i\right)i=1,2,\cdots ,k$. The equation of the fitted space line is:

$$\frac{x-x_0}{l}=\frac{y-y_0}{m}=\frac{z-z_0}{n}$$

(1)

where $\left(l,m,n\right)$ is the direction vector of the space straight line, and $\left(x_0,y_0,z_0\right)$ is the coordinate of a point on the space straight line. Equation (1) is converted into a spatial straight line projective equation:

$$\left\{\begin{array}{l}x=\frac{l}{n}(z-z_0)+x_0=a_{l}z+b_l\\ y=\frac{m}{n}(z-z_0)+y_0=c_{l}z+d_l\end{array}\right.$$

(2)

where $a_l=\frac{l}{n}$, $b_l=x_0-\frac{l}{n}{z}_0$, $c_l=\frac{m}{n}$, $d_l=y_0-\frac{m}{n}{z}_0$.

The observation equations are constructed to satisfy:

$$\left\{\begin{array}{c}{Q}_x={\displaystyle \sum _{i=1}^k{\left[x_i-\left(a_l{z}_i+b_l\right)\right]}^2=\min }\\ Q_y={\displaystyle \sum _{i=1}^k{\left[y_i-\left(c_l{z}_i+d_l\right)\right]}^2=\min }\end{array}\right.$$

(3)

Based on the data of $k$ 3D coordinates, the values of $a_l,b_l,c_l,d_l$ are estimated by solving the equations with the LS, and the gravity direction vector $g_{LTCS}$ of the free-falling trajectory in LTCS is obtained as $\left(\begin{array}{ccc}{a}_l& c_l& 1\end{array}\right)$.

2.1.2. Spatial Linear Fitting Using TLS

The spatial straight line projective Equation (2) is transformed into matrix form [30]:

$$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{cccc}z& 1& 0& 0\\ 0& 0& z& 1\end{array}\right]\left[\begin{array}{c}{a}_l\\ b_l\\ c_l\\ d_{}\end{array}\right]$$

(4)

According to Equation (4), the error equation is constructed:

$$V_l=\left[\begin{array}{cccc}z& 1& 0& 0\\ 0& 0& z& 1\end{array}\right]\left[\begin{array}{c}{\widehat{a}}_l\\ {\widehat{b}}_l\\ {\widehat{c}}_l\\ {\widehat{d}}_l\end{array}\right]-\left[\begin{array}{c}x\\ y\end{array}\right]$$

(5)

Let ${{\rm B}}_l=\left[\begin{array}{cccc}z& 1& 0& 0\\ 0& 0& z& 1\end{array}\right]$, $L_l=\left[\begin{array}{c}x\\ y\end{array}\right]$,${\widehat{X}}_l={\left[\begin{array}{cccc}{\widehat{a}}_l& {\widehat{b}}_l& {\widehat{c}}_l& {\widehat{d}}_l\end{array}\right]}^T$, and then Equation (5) is simplified to:

$$V_l={{\rm B}}_l{\widehat{X}}_l-L_l$$

(6)

Based on the 3D observations $L_i$ of the spatial line, the error Equation (6) is solved using the TLS iterative method to meet the adjustment criteria:

$${\displaystyle \sum _{i=1}^k{\left({\widehat{L}}_i-L_i\right)}^2+{\displaystyle \sum _{j=1,i=1}^{j=t,i=k}{\left({\widehat{B}}_l^{ij}-B_l^{ij}\right)}^2}}=\min$$

(7)

By substituting Equation (6) and taking derivatives of each element, the iterative equation is obtained:

$$\left\{\begin{array}{c}{\widehat{B}}_l^T{\widehat{X}}_l={\widehat{B}}_l^T{L}_l\\ N_b{\widehat{B}}_l^T=B_l^T+{\widehat{X}}_l{L}_l^T\end{array}\right.$$

(8)

where $N_b=E+{\widehat{X}}_l{\widehat{X}}_l^T$. The algorithm flow of the TLS iterative method is shown in Figure 3.

2.1.3. Comparative Analysis of LS and TLS

To illustrate the accuracy of the two methods for extracting gravity direction vectors, we carry out simulation experiments. Assuming that the point precision of the laser tracker is 10 μm, the falling distance of the free-fall target is about 1 m, and the laser tracking measurement frequency can reach 2000 Hz. The theoretical accuracy of the two methods for extracting gravity direction vectors is calculated by 1000 free-falling simulation experiments. The average of 30 observations was recorded as one observation, and the results are shown in Figure 4 and

Table 1. Statistics of gravity direction vectors of LS and TLS (″).

	MAX	MIN	MEAN	STD	RMS
LS	0.192	0.102	0.146	0.014	0.146
TLS	0.192	0.102	0.146	0.014	0.146

.

From Figure 4 and Table 1, it can be seen that the STD and RMS of LS and TLS are the same, both better than 0.15, indicating that the accuracy of the two methods is consistent. However, since the TLS involves an iterative calculation, which is more complicated and time-consuming compared to LS, the LS is chosen for extracting gravity direction vector in this paper.

2.2. Space Direction Vector Transformation

In the proposed DOV measurement method, it is necessary to convert the gravity direction vector in LTCS to GCCS, but this transformation involves the rotation of the large rotation angle, and the traditional Bursa model is no longer applicable. Based on the nonlinear 13-parameter transformation model with a large rotation angle [31,32], this paper proposes a direction transformation method using baseline information to realise the transformation of the gravity direction vector in LTCS to GCCS.

Several GNSS antennas are fixed onto the measuring device, and the coordinates of the antenna phase centre in LTCS are measured using the laser tracker, and the distance between each antenna is calculated to obtain the baseline information. Because the measurement precision of the laser tracker reaches the sub-micrometre level, the sub-micrometre precision distance constraint is added to GNSS baseline solution to obtain high-precision GNSS baselines. Based on several baselines, the direction transformation parameters are solved as shown below.

Based on the nonlinear 13-parameter transformation model, the relationship between the common point $A_i$ in GCCS and LTCS can be expressed as:

$$\left[\begin{array}{l}{X}_{A_i}\\ Y_{A_i}\\ Z_{A_i}\end{array}\right]=\psi R\left[\begin{array}{l}{x}_{A_i}\\ y_{A_i}\\ z_{A_i}\end{array}\right]+\left[\begin{array}{l}Tx\\ Ty\\ Tz\end{array}\right]i=1,2,\cdots ,n$$

(9)

where $\left(x_{A_i},y_{A_i},z_{A_i}\right)$ is the coordinate of the common point $A_i$ in LTCS, $\left(X_{A_i},Y_{A_i},Z_{A_i}\right)$ is the coordinate of the common point $A_i$ in GCCS, $R=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$, $\left(a_1,a_2,a_3\right),\left(b_1,b_2,b_3\right),\left(c_1,c_2,c_3\right)$ are the direction cosines of the x, y, and z axes in GCCS, $\psi$ is the scale parameter, and $\left(Tx,Ty,Tz\right)$ are the translation parameters.

Since the rotation matrix $R$ is an orthogonal matrix and the corresponding coordinate transformation is an orthogonal transformation, the conditional equation is obtained:

$$\left\{\begin{array}{l}{a}_1{}^2+a_2{}^2+a_3{}^2=1\\ b_1{}^2+b_2{}^2+b_3{}^2=1\\ c_1{}^2+c_2{}^2+c_3{}^2=1\\ a_1{a}_2+b_1{b}_2+c_1{c}_2=0\\ a_1{a}_3+b_1{b}_3+c_1{c}_3=0\\ a_2{a}_3+b_2{b}_3+c_2{c}_3=0\end{array}\right.$$

(10)

From Equation (10), it can be seen that there are only three independent parameters in the rotation matrix $R$, taking $a_2,a_3,b_3$ as the independent parameters. The remaining six parameters are expressed as:

$$\left\{\begin{array}{l}{a}_1=\sqrt{1-a_2{}^2-a_3{}^2}\\ c_3=\sqrt{1-a_3{}^2-b_3{}^2}\\ b_1=\frac{-a_1{a}_3{b}_3-a_2{c}_3}{1-a_3{}^2}\\ b_2=\sqrt{1-b_1{}^2-b_3{}^2}\\ c_1=a_2{b}_3-a_3{b}_2\\ c_2=a_3{b}_1-a_1{b}_3\end{array}\right.$$

(11)

According to Equation (9), the common points $A_i$ and $A_j$ respectively construct equations. The two equations are differenced to eliminate the translation parameters. The scale parameters will be added to each baseline separately:

$$\left[\begin{array}{l}\Delta X_{A_{ij}}\\ \Delta Y_{A_{ij}}\\ \Delta Z_{A_{ij}}\end{array}\right]={\psi }_{ij}R\left[\begin{array}{l}\Delta x_{A_{ij}}\\ \Delta y_{A_{ij}}\\ \Delta z_{A_{ij}}\end{array}\right]$$

(12)

where ${\psi }_{ij}$ is the scale parameter of the baseline $A_{ij}$, $\left(\begin{array}{ccc}\Delta X_{Aij}& \Delta Y_{A_{ij}}& \Delta Z_{A_{ij}}\end{array}\right)$ are the GNSS baselines of common points $A_i$ and $A_j$, and $\left(\begin{array}{ccc}\Delta x_{A_{ij}}& \Delta y_{Aij}& \Delta z_{Aij}\end{array}\right)$ are the common point baselines in LTCS.

According to more than three common baselines, the rotation matrix $R$ can be obtained by solving Equation (12) using the LS method. However, the rotation matrix $R$ only has three independent parameters, and the remaining six parameters are its nonlinear functions. It is very complicated to directly solve Equation (12), so the first-order expansion of Taylor series is used to solve it.

Assume that the unknowns are nine direction cosine parameters, $N$ scale parameters, and $N$ is the number of common baselines. Then, Equation (12) is expanded by Taylor series to obtain:

$$\begin{array}{c}\left[\begin{array}{c}\Delta X_{A_{ij}}\hfill \\ \Delta Y_{A_{ij}}\hfill \\ \Delta Z_{A_{ij}}\hfill \end{array}\right]={\psi }_{{}_{ij}}^0\left[\begin{array}{ccc}{a}_1^0& a_2^0& a_3^0\\ b_1^0& b_2^0& b_3^0\\ c_1^0& c_2^0& c_3^0\end{array}\right]{\left[\begin{array}{c}\Delta x_{A_{ij}}\hfill \\ \Delta y_{A_{ij}}\hfill \\ \Delta z_{A_{ij}}\hfill \end{array}\right]}_i+\left[\begin{array}{c}{a}_1^0\Delta x_{A_{ij}}+a_2^0\Delta y_{A_{ij}}+a_3^0\Delta z_{A_{ij}}\hfill \\ b_1^0\Delta x_{A_{ij}}+b_2^0\Delta y_{A_{ij}}+b_3^0\Delta z_{A_{ij}}\hfill \\ c_1^0\Delta x_{A_{ij}}+c_2^0\Delta y_{A_{ij}}+c_3^0\Delta z_{A_{ij}}\hfill \end{array}\right]\mathrm{d}{\psi }_{ij}+\hfill \\ \left[\begin{array}{ccccccccc}{\psi }_{{}_{ij}}^0\Delta x_{A_{ij}}& {\psi }_{{}_{ij}}^0\Delta y_{A_{ij}}& {\psi }_{{}_{ij}}^0\Delta z_{A_{ij}}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\psi }_{{}_{ij}}^0\Delta x_{A_{ij}}& {\psi }_{{}_{ij}}^0\Delta y_{A_{ij}}& {\psi }_{{}_{ij}}^0\Delta z_{A_{ij}}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\psi }_{{}_{ij}}^0\Delta x_{A_{ij}}& {\psi }_{{}_{ij}}^0\Delta y_{A_{ij}}& {\psi }_{{}_{ij}}^0\Delta z_{A_{ij}}\end{array}\right]\hfill \\ {\left[\begin{array}{ccccccccc}\mathrm{d}{a}_1& \mathrm{d}{a}_2& \mathrm{d}{a}_3& \mathrm{d}{b}_1& \mathrm{d}{b}_2& \mathrm{d}{b}_3& \mathrm{d}{c}_1& \mathrm{d}{c}_2& \mathrm{d}{c}_3\end{array}\right]}^{\mathrm{T}}\hfill \end{array}$$

(13)

where the variables with superscript 0 represent the approximate value of the variable, and $\mathrm{d}{\psi }_{ij}$,$\mathrm{d}{a}_1$,$\mathrm{d}{a}_2$,$\mathrm{d}{a}_3$,$\mathrm{d}{b}_1$,$\mathrm{d}{b}_2$,$\mathrm{d}{b}_3$,$\mathrm{d}{c}_1$,$\mathrm{d}{c}_2$,$\mathrm{d}{c}_3$ are the corrections. Equation (13) is converted into the form of the error equation to obtain:

$$V_{dt}=A_{dt}{X}_{dt}-L_{dt}$$

(14)

where $V_{dt}^{ij}={\left[\begin{array}{ccc}{V}_{\Delta X_{A_{ij}}}& V_{\Delta Y_{A_{ij}}}& V_{\Delta Z_{A_{ij}}}\end{array}\right]}^{\mathrm{T}}$ denotes the corrections of each GNSS baseline $A_{ij}$ in the X, Y, and Z directions, respectively, and $X_{dt}={\left[\begin{array}{llllllllll}\mathrm{d}{\psi }_{ij}\hfill & \mathrm{d}{a}_1\hfill & \mathrm{d}{a}_2\hfill & \mathrm{d}{a}_3\hfill & \mathrm{d}{b}_1\hfill & \mathrm{d}{b}_2\hfill & \mathrm{d}{b}_3\hfill & \mathrm{d}{c}_1\hfill & \mathrm{d}{c}_2\hfill & \mathrm{d}{c}_3\hfill \end{array}\right]}^{\mathrm{T}}$.

According to Equation (11), the conditional equation is listed to obtain:

$$A_{dt}^{\prime }{X}_{dt}+W_{dt}=0$$

(15)

where $X_{dt}$ means the same as above, and $A_{dt}^{\prime }$ and $W_{dt}$ are:

$$A_{dt}^{\prime }=\left[\begin{array}{cccccccccc}0\hfill & 2a_1^0& 2a_1^0& 2a_1^0& 0& 0& 0& 0& 0& 0\\ 0\hfill & 0& 0& 0& 2b_1^0& 2b_2^0& 2b_3^0& 0& 0& 0\\ 0\hfill & 0& 0& 0& 0& 0& 0& 2c_1^0& 2c_2^0& 2c_3^0\\ 0\hfill & a_2^0& a_1^0& 0& b_2^0& b_1^0& 0& c_2^0& c_1^0& 0\\ 0\hfill & a_3^0& 0& a_1^0& b_3^0& 0& b_1^0& c_3^0& 0& c_1^0\\ 0\hfill & 0& a_3^0& a_2^0& 0& b_3^0& b_2^0& 0& c_3^0& c_2^0\end{array}\right];$$

$$W_{dt}=\left[\begin{array}{l}{a}_1^0+a_2^0+a_3^0-1\\ b_1^0+b_2^0+b_3^0-1\\ z_1^0+z_2^0+z_3^0-1\\ a_1^0{a}_2^0+b_1^0{b}_2^0+c_1^0{c}_2^0\\ a_1^0{a}_3^0+b_1^0{b}_3^0+c_1^0{c}_3^0\\ a_2^0{a}_3^0+b_2^0{b}_3^0+c_2^0{c}_3^0\end{array}\right].$$

To improve the calculation efficiency, the conditional equation is transformed into a pseudo-observation equation, and then a new error equation is formed with Equation (14):

$$\left[\begin{array}{c}{V}_{dt}\\ V_{dt}^{\prime }\end{array}\right]=\left[\begin{array}{c}{A}_{dt}\\ A_{dt}^{\prime }\end{array}\right]X_{dt}+\left[\begin{array}{c}{L}_{dt}\\ W_{dt}\end{array}\right]$$

(16)

The parameters $X_{dt}$ are obtained by solving Equation (11) using the indirect adjustment.

The iterative computation is performed to obtain the rotation matrix $R$, and the steps are as follows:

(i)

The approximate value is generally desirable:

$${\psi }_{ij}=1,R=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].$$

(17)

(ii)

The error equation is composed according to Equation (11), and if there are $N$ baselines, $3N+6$ error equations can be composed.

(iii)

The corrections of $N+9$ unknowns are solved by $3N+6$ equations.

(iv)

The latest values of unknowns are calculated.

(v)

According to the corrections, judge whether the convergence requirement is satisfied.

If not, repeat steps (ii) to (v) until the convergence is satisfied.

(vi)

According to the rotation matrix $R_T$ obtained iteratively, the gravity direction vector $g_{LTCS}$ in LTCS is converted to GCCS:

$$g_{GCCS}=R_T{g}_{LTCS}$$

(18)

where $g_{LTCS}$ is the gravity direction vector in GCCS.

2.3. Calculation of Astrogeodetic DOV

The astrogeodetic DOV, also known as the relative DOV, has relative significance because the normals of the measured points under different reference ellipsoids are different, and the astrogeodetic DOV of the point is also different. Users can obtain the DOV of the corresponding ellipsoid through different ellipsoid transformations. The WGS1984 ellipsoid is selected as an example to illustrate the principle of the proposed method in the paper.

The semi-major axis of the reference ellipsoid corresponding to GCCS is $\mathrm{a}$. The earth’s ellipticity is $f$, and the semi-minor axis of the ellipsoid is $\mathrm{b}=\mathrm{a}-\mathrm{a}*f$. The reference ellipsoid can be expressed as:

$$\frac{X^2}{{\mathrm{a}}^2}+\frac{Y^2}{{\mathrm{a}}^2}+\frac{Z^2}{{\mathrm{b}}^2}=1$$

(19)

Let $F\left(X,Y,Z\right)=\frac{X^2}{{\mathrm{a}}^2}+\frac{Y^2}{{\mathrm{a}}^2}+\frac{Z^2}{{\mathrm{b}}^2}-1$, and for X, Y, and Z, the first order derivatives are:

$$\left\{\begin{array}{c}{F}_X\left(X,Y,Z\right)=\frac{2}{{\mathrm{a}}^2}\cdot X\\ F_Y\left(X,Y,Z\right)=\frac{2}{{\mathrm{a}}^2}\cdot Y\\ F_Z\left(X,Y,Z\right)=\frac{2}{{\mathrm{b}}^2}\cdot Z\end{array}\right.$$

(20)

Assuming that the geodetic cartesian coordinates of the measuring point P measured by GNSS are $\left(X_0,Y_0,Z_0\right)$, the ellipsoid normal vector through point P is:

$$p=\left(\frac{2}{{\mathrm{a}}^2}\cdot X_0,\frac{2}{{\mathrm{a}}^2}\cdot Y_0,\frac{2}{{\mathrm{b}}^2}\cdot Z_0\right)$$

(21)

In GCCS, the meridian plane through point P can be determined by the space vectors $(0,0,1)$ and $(X_0,Y_0,Z_0)$. Assuming that the normal direction vector of the meridian plane is $\left(X_{me},Y_{me},Z_{me}\right)$, there are:

$$\left\{\begin{array}{l}{Z}_{me}=0\\ X_0{X}_{me}+Y_0{Y}_{me}+Z_0{Z}_{me}=0\end{array}\right.$$

(22)

Let $X_{me}=1$, and the normal direction vector of the meridian plane through point P is $p_{me}\left(1,-\frac{X_0}{Y_0},0\right)$.

The prime vertical plane through point P can be determined by the space vectors $p\left(\frac{2}{a^2}\cdot X_0,\frac{2}{a^2}\cdot Y_0,\frac{2}{b^2}\cdot Z_0\right)$ and $p_{me}\left(1,-\frac{X_0}{Y_0},0\right)$. Assuming that the normal direction vector of the prime vertical plane is $\left(X_{pr},Y_{pr},Z_{pr}\right)$, there are:

$$\left\{\begin{array}{l}\frac{2}{a^2}{X}_0{X}_{pr}+\frac{2}{a^2}{Y}_0{Y}_{pr}+\frac{2}{b^2}{Z}_0{Z}_{pr}=0\\ X_{pr}-\frac{X_0}{Y_0}{Y}_{pr}=0\end{array}\right.$$

(23)

Let $X_{pr}=1$, and the normal direction vector of the prime vertical plane through point P is $p_{pr}\left(1,\frac{Y_0}{X_0},-\frac{b^2}{a^2}\cdot \frac{X_0{}^2+Y_0{}^2}{X_0{Z}_0}\right)$.

According to the gravity direction vector $g_{GCCS}$ and ellipsoid normal vector $p$, the astrogeodetic DOV $u$ is calculated as:

$$u=\arccos \frac{g_{GCCS}\cdot p}{\left|g_{GCCS}\right|\left|p\right|}$$

(24)

The gravity direction vector $g_{GCCS}$ is projected onto the meridian and prime vertical planes, respectively:

$$\left\{\begin{array}{l}{g}_{me}=g_{GCCS}-\frac{p_{me}}{\left|p_{me}\right|}\cdot g_{GCCS}\cdot \frac{p_{me}}{\left|p_{me}\right|}\\ g_{pr}=g_{GCCS}-\frac{p_{pr}}{\left|p_{pr}\right|}\cdot g_{GCCS}\cdot \frac{p_{pr}}{\left|p_{pr}\right|}\end{array}\right.$$

(25)

where $g_{me}$ is the meridian component of $g_{GCCS}$, and $g_{pr}$ is the prime vertical component of $g_{GCCS}$.

The meridian and prime vertical components of astrogeodetic DOV are calculated as:

$$\left\{\begin{array}{l}\xi =\arccos \frac{g_{me}\cdot p}{\left|g_{me}\right|\left|p\right|}\\ \eta =\arccos \frac{g_{pr}\cdot p}{\left|g_{pr}\right|\left|p\right|}\end{array}\right.$$

(26)

where $\xi$ is the meridian component of astrogeodetic DOV, and $\eta$ is the prime vertical component of astrogeodetic DOV.

According to the meridian and prime vertical components of astrogeodetic DOV at the measuring point P, and the geodetic coordinate $\left(B,L\right)$, the astronomical longitude and latitude $\left(\phi ,\lambda \right)$ are calculated using the following relationship:

$$\left\{\begin{array}{l}\phi =B+\xi \\ \lambda =L+\eta \sec \phi \end{array}\right.$$

(27)

3. Simulation Study

3.1. Theoretical Precision of DOV

The main error sources of the proposed method are laser tracker measurement error and GNSS measurement error. The influence of laser tracker error on the DOV measurement is reflected in two aspects: one is that when tracking the free-falling target, the measurement errors will make the measured gravity direction vector different from the theoretical vector. Second, the errors of the laser tracker will affect the common point coordinates during the direction transformation between the GCCS and LTCS, and then affect the direction transformation parameters. Similarly, the GNSS baselines are used to solve the direction transformation parameters between the GCCS and LTCS, and its solution precision determines the accuracy of the rotation matrix. The GNSS coordinate precision of the measurement point P also affects the ellipsoidal normal vector, but its effect on the ellipsoidal normal vector is negligible due to the large lengths of the semi-major and semi-minor axes of the earth.

The simulation experiments are carried out to analyse the theoretical precision of the proposed measurement method. Assuming that the point precision of the laser tracker is 10 μm, the falling distance of the free-fall target is about 1 m, and the laser tracking measurement frequency can reach 2000 Hz, so the observation time of each free-fall measurement is about 0.45 s, and the sampling data reaches about 900. Before observing the free fall, the fixed GNSS antennas were measured 1000 times to determine the coordinates in LTCS. Suppose the GNSS positioning precision is 1 mm, and there are four GNSS devices in total. Due to the addition of the distance constraints of sub-micrometre precision during the GNSS baseline solving, the calculated GNSS baseline accuracy is set to 0.1 mm. The reference geodetic cartesian coordinates of measuring point P are: X = −2,148,744.4656 m, Y = 4,426,641.1849 m, and Z = 4,044,656.0516 m. The ellipsoid WGS1984 is selected, the ellipsoid semi-major axes a = 6,378,137 m, and the earth’s ellipticity $f$ = 1/298.257223563.

Firstly, a 3D coordinate series of the falling free-fall target is simulated, and the laser tracker error is added to the theoretical coordinates to generate the simulated observations in LTCS. The method using the fixed GNSS antennas to obtain the direction transformation parameters of GCCS and LTCS is as follows: the theoretical values of the fixed GNSS antennas in LTCS are known in advance, and the laser tracker error is added to the theoretical coordinates to generate simulation observations of the fixed GNSS antennas in LTCS. The theoretical coordinates of the fixed GNSS antennas in LTCS and the theoretical transformation parameters between the two coordinate systems are used to obtain the theoretical coordinates of the fixed GNSS antennas in GCCS, and then the theoretical coordinates of each GNSS antenna are differenced to obtain the theoretical values of the baseline vectors. According to the solving precision of the GNSS baseline, the error is added to the GNSS baseline vector to obtain its simulated observations. Based on the simulated observations of the fixed GNSS antennas in LTCS and the simulated observation baseline vectors in GCCS, the direction transformation parameters between the two coordinate systems affected by the errors of the laser tracker and GNSS can be obtained using the space direction vector transformation method in Section 2.3. The direction transformation parameters with errors between GCCS and LTCS are used to calculate the gravity direction vector containing errors in GCCS. Combined with the ellipsoid normal vector determined by the geodetic coordinates of the measuring point P, the astrogeodetic DOV with errors is calculated. Finally, the theoretical precision of DOV can be obtained by calculating the difference between the error-containing DOV and the theoretical DOV.

Four GNSS antennas are evenly distributed in a circle with the radius r centred on the laser tracker. To consider the size of the measuring device, we separately discussed the precision of the DOV for different radius r. The average of 30 observations was recorded as one observation, and the above experiments were repeated 1000 times. The results are shown in Figure 5 and

Table 2. Precision statistics of DOV with different radius r (″).

	Radius r	2 m	2.5 m	3 m	3.5 m
u	MAX	0.91	0.78	0.63	0.72
	MIN	−0.85	−0.75	−0.72	−0.48
	MEAN	0.00	−0.01	−0.01	0.01
	STD	0.28	0.24	0.20	0.17
	RMS	0.28	0.24	0.20	0.17
ξ	MAX	1.13	0.91	0.69	0.76
	MIN	−0.99	−0.85	−0.84	−0.65
	MEAN	−0.01	−0.01	−0.01	0.02
	STD	0.32	0.28	0.23	0.20
	RMS	0.32	0.28	0.23	0.20
η	MAX	1.12	0.88	0.83	0.54
	MIN	−1.31	−0.93	−0.65	−0.55
	MEAN	0.00	0.00	0.00	0.00
	STD	0.33	0.25	0.22	0.19
	RMS	0.33	0.25	0.22	0.19

.

It can be seen from Figure 5 and Table 2 that the proposed measurement method has good stability, the error fluctuates between ±1.2″, and the average value approaches 0. With the increase of the radius r of the circle formed by four GNSS antennas, the measurement precision of the DOV becomes higher, which basically meets the requirements of the first-class astronomical precision (0.3″) in China. When the radius r is 3.5 m, the precision can reach 0.2″.

3.2. Influence of Laser Tracking Measurement Error on DOV

To discuss the influence of the laser tracking measurement error on the DOV, it is assumed that the GNSS baseline solution is error-free, and only the laser tracking measurement error is added to the simulated theoretical values, and the remaining conditions are consistent with Section 3.1. The above experiments were repeated 1000 times, and the results are shown in Figure 6 and

Table 3. Influence statistics of laser tracking error on DOV (″).

	Radius r	2 m	2.5 m	3 m	3.5 m
u	MAX	0.07	0.07	0.07	0.06
	MIN	−0.07	−0.07	−0.07	−0.06
	MEAN	0.00	0.00	0.00	0.00
	STD	0.02	0.02	0.02	0.02
	RMS	0.02	0.02	0.02	0.02
ξ	MAX	0.07	0.07	0.07	0.07
	MIN	−0.06	−0.07	−0.07	−0.07
	MEAN	0.00	0.00	0.00	0.00
	STD	0.02	0.02	0.02	0.02
	RMS	0.02	0.02	0.02	0.02
η	MAX	0.07	0.07	0.07	0.06
	MIN	−0.08	−0.07	−0.07	−0.07
	MEAN	0.00	0.00	0.00	0.00
	STD	0.02	0.02	0.02	0.02
	RMS	0.02	0.02	0.02	0.02

.

It can be seen from Figure 6 and Table 3 that the influence of the laser tracking measurement error on the DOV is ±0.07″, the STD of the DOV is 0.02″, and the effect remains consistent for each direction and does not change as the circle radius r increases. In general, the influence of laser tracking measurement on the DOV accounts for about 2% of the total level, which is relatively stable.

3.3. Influence of GNSS Error on DOV

To discuss the influence of the GNSS error on DOV, it is assumed that the laser tracking measurement error is 0, and only the GNSS error is added to the simulated theoretical value of the GNSS baseline, and the remaining conditions are consistent with Section 3.1. The above experiments were repeated 1000 times, and the results are shown in Figure 7 and

Table 4. Influence statistics of GNSS error on DOV (″).

	Radius r	2 m	2.5 m	3 m	3.5 m
u	MAX	0.81	0.73	0.58	0.50
	MIN	−0.83	−0.69	−0.57	−0.57
	MEAN	0.01	0.00	0.00	0.00
	STD	0.29	0.25	0.20	0.17
	RMS	0.29	0.25	0.20	0.17
ξ	MAX	1.00	0.86	0.68	0.59
	MIN	−1.04	−0.86	−0.63	−0.68
	MEAN	0.00	0.00	0.00	0.00
	STD	0.33	0.28	0.23	0.20
	RMS	0.33	0.28	0.23	0.20
η	MAX	1.05	0.94	0.61	0.53
	MIN	−0.98	−0.94	−0.74	−0.48
	MEAN	0.01	0.00	0.00	0.01
	STD	0.32	0.25	0.21	0.18
	RMS	0.32	0.25	0.21	0.18

.

As can be seen from Figure 7 and Table 4, the GNSS baseline solution error is the decisive factor affecting the DOV, and the single influence of GNSS is basically consistent with the simulation precision of the proposed method. With the increase of the circle radius r formed by the four GNSS antennas, the measurement precision of the DOV increases gradually, which is due to the increase of the baseline length, resulting in the decrease of the relative influence of GNSS error. Overall, GNSS error is the key factor for the all-weather DOV measurement method proposed in the paper.

4. Conclusions

A novel all-weather method to determine the DOV by combining 3D laser tracking free-fall and multi-GNSS baselines was proposed. The laser tracker was used to dynamically observe the three-dimensional coordinate series of the free-falling target at the measuring point, and the gravity direction vector in the laser tracker coordinate system was obtained by the least squares fitting coordinate series. In this paper, a direction transformation method using baseline information was proposed to convert the gravity direction vector in the laser tracker coordinate system to the geodetic cartesian coordinate system. The ellipsoid normal vector at the measuring point can be calculated by the derived formula from the geodetic coordinates obtained by GNSS measurement, and then the astrogeodetic DOV of the measuring point can be estimated. The DOV was projected to the meridian and prime vertical planes to obtain the meridian and prime vertical components of the DOV, respectively. The astronomical latitude and longitude can be calculated according to meridian and prime vertical components of DOV and geodetic coordinates. Through theoretical analysis, the method has good feasibility.

The main error sources of the proposed DOV measurement method are laser tracking and GNSS measurement errors. The control variable method was adopted to discuss its influence on the DOV. The results showed that the contribution of the laser tracking measurement error on the DOV was only 2%, and remained stable. The GNSS error was the key factor of the proposed all-weather DOV measurement method. Overall, the proposed method can simultaneously determine the astrogeodetic DOV and the astronomical longitude and latitude at the measuring point. Compared with the traditional methods of astronomical geodesy and astronomical gravimetry, the proposed method is hardly affected by the climate and environment, it can be observed in all weathers, and has a high theoretical precision of 0.2″.