A Conical Model Approach for Invariant Points of Very Long Baseline Interferometry and Satellite Laser Ranging

Abstract

A new realization of the international terrestrial reference frame, a combination of four different space geodetic techniques, was released in 2022. Each geodetic solution should be combined carefully based on the local tie information at the co-located site. Although many approaches have been successfully applied to connect different geodetic sensors, to date, there has been no unified mathematical representation for the target motions. Herein, a unified conical model was developed to estimate the invariant points of geodetic sensors using a more robust and consistent approach. It modeled the motion of targets, in either the horizontal or vertical axis, as cones; thus, homogeneous modeling was implemented. In addition to its simplicity, the model simultaneously estimated the tilting of the vertical axis and horizontal offset. The mathematical relationship and normality of the normal vector were modeled as a Gauss–Markov model with fixed constraints. The pre-computed initial coordinates of the pillars and targets were adjusted simultaneously to calculate the correlation information for the local tie vector. The complete model was successfully applied to the co-located site, which was transformed into a global reference frame via the Helmert transformation based on the global navigation satellite system campaign. The results showed that the proposed method is more efficient in terms of the number of parameters for invariant points of geodetic sensors (only 13% compared to the 3D circle fitting type conventional approach). In addition, the reliability of the estimated solution can be increased by avoiding an ill-conditioned linear system through the conical model.

1. Introduction

The International Earth Rotation Service (IERS) plays a crucial role in establishing and maintaining the international terrestrial reference frame (ITRF) by analyzing data from four space geodetic systems as follows: global navigation satellite system (GNSS), satellite laser ranging (SLR), very long baseline interferometry (VLBI), and Doppler orbitography and radiopositioning integrated by satellites [1,2,3,4,5,6,7,8]. The IERS network is defined using space geodetic instruments installed on the sites, and many of the stations have more than two geodetic instruments located very close to each other [9]. The “co-located site” is the station where more than one instrument is in operation simultaneously. It is recommended to link the reference points of instruments using a local tie vector that represents a 3D baseline vector between the reference points of two space geodesy instruments [10,11]. Therefore, it is necessary to estimate a local tie vector with high accuracy to establish and maintain the geodetic reference frame, i.e., ITRF.

The reference point (or invariant point, IVP) of geodetic instruments, such as the VLBI telescope and SLR transmitter, is defined as the projection point of the secondary (elevation, EL) axis onto the primary (azimuth, AZ) axis of rotation [12]. Generally, the IVP of an instrument is located inside it and is thus not accessible. Therefore, the indirect method is frequently used to determine the location of the IVP [13,14,15] and involves a terrestrial ground survey by observing the reflective targets attached to the surface of the instrument. Subsequently, the trace of each target created by the rotational motions of the primary and secondary axes is modeled as a 3D circle in space [10,11,13,14,16,17,18,19,20]. Various types of 3D circle fitting algorithms (called “conventional” methods) have been proposed to compute the location of the IVP using terrestrial ground survey techniques. However, only a partial arc, instead of a full circle, can be obtained from the rotational sequences of the elevation axis because of the limited rotational motion of the elevation axis. This leads to an ill-conditioned linear system; consequently, the quality of the estimated IVP is degraded. To overcome this, modeling an incomplete arc as a part of the cone, i.e., its base, has been attempted. This approach significantly reduces the number of estimated parameters and increases the estimability of the linear system solved using the least-squares method [21]. However, the cone-shape model can only be applied to arcs generated from the rotational motion about the elevation axis and has a limitation in the estimation of the possible offset between the azimuth and elevation axes and the tilt angle between the azimuth axis and zenith direction.

In this study, the conical model is applied not only to the arcs from the rotational motion about the elevation axis but also to circles from the rotational motion about the azimuth axis. The conical model approach has the following advantages: (1) immaculate consistency in setting up the mathematical functions for both arcs from the rotational motions and (2) flexibility in the modeling, although tilting of the azimuth axis may exist. As mentioned earlier, IVP is defined as the point at which the elevation axis is projected onto the azimuth axis when an offset exists between the two axes. Hence, the mathematical model was fully extended to describe the possible offset between the azimuth and elevation axes. Once the complete mathematical model is constructed to connect the parameters with measurements from terrestrial ground surveying, adjustment computation is performed to obtain the least-squares solution. The new conical model provides consistency and efficiency in terms of geometric relationship for the IVP and numerical stability of the parameters with the corresponding covariance information.

In the following sections, detailed descriptions of the measurements from terrestrial ground surveying techniques are introduced. Then, a mathematical model to accommodate the tilting of the azimuth axis and offset between the axes and an adjustment model for the least-squares solution is proposed. Subsequently, the numerical test results and conclusions are presented.

2. Local Tie Strategy

Recently, a new realization of the international terrestrial reference system (known as ITRF2020) was announced by combining the individual solutions of space geodetic techniques [2]. To combine different solutions, the relative positions between the reference points should be represented together at the co-located sites. Because the reference points cannot be accessed directly, a special procedure for estimating the IVPs and the corresponding correlation information is required [22].

One of the co-located sites is operational in Korea (Sejong station), of which VLBI and GNSS are maintained by the National Geographic Information Institute, while SLR is operated by the Korea Astronomy and Space Science Institute. The diameter of the VLBI antenna is 22 m, and the GNSS station (4-char site name SEJN) is one of the international GNSS service (IGS) continuously operating reference stations (CORS). As shown in Figure 1, the two sensors are approximately 200 m apart and surrounded by pillars (five pillars for VLBI and three for SLR). The ground surveying to connect these sensors is accomplished by measuring the slant distances, the horizontal/vertical angles, and the height differences by spirit leveling; however, only two VLBI pillars are involved because of the obstruction by the nearby trees. Thus, the geometric condition of the connecting two areas is weak in rotation parameters.

The observations can be divided into two surveying operations as follows: (1) conventional surveying in local frames, such as slant distances, horizontal and vertical angles, and height differences, and (2) GNSS campaign in the global reference frame. For local surveying, a total station was installed on the pillars, except at the GNSS station, with the origin defined at VP01 (VP02 is defined on the $\mathit{x}$-axis or east direction). A 7-day GNSS campaign for all nine pillars, including the GNSS station, was conducted to calculate the coordinates of the global reference frame. The collected GNSS data were processed simultaneously with the IGS core stations to align with the global reference frame.

The local tie vectors (GNSS to IVPs), along with their covariance matrices, are estimated in the local reference frame; thus, they should be transformed into the global frame. The Helmert transformation is generally applied to estimate the parameters (translation, rotation, and scale) between the reference frames, assuming small rotation angles [23].

The developed mathematical model is based on two types of targets, either horizontal or vertical, although several of the targets are involved in both motions, depending on the attached location. Both VLBI and SLR have their own azimuthal systems that are not identical to the local geodetic reference frame. It should be mentioned here that we aim at the azimuth-elevation (AZ-EL) type, which is the most common mount among several types of VLBI antennas [24].

Note that the rotational axes for both motions are independent of each other. More specifically, the horizontal rotation ranges from $0^{\circ }$ to ${360}^{\circ }$ in (its own) azimuthal system, with a zenith angle of $0^{\circ }$ at all times, whereas the vertical motion covers only ${90}^{\circ }$ (i.e., a quarter circle), as shown in Figure 2. Therefore, the geometric condition should be reflected by considering the characteristics of the motion, as described in the following section. Detailed target location and ground survey information can be found in [21].

3. The Conical Model

The sphero-conical model was successfully applied to estimate the VLBI IVP [21]. However, the horizontal and vertical motions, i.e., the circles in 3D space (horizontal targets) and the cones with the apex at the IVP (vertical targets), respectively, were modeled differently. Therefore, this model is insufficient for universal applications and has some limitations in that the $\mathit{z}$-component of the IVP is not involved in the horizontal motion of the targets. Furthermore, if the vertical axis is tilted by chance, the angle cannot be readily calculated using the sphero-conical model.

Instead of separately modeling both motions, the model was modified to apply a homogeneous geometrical relationship (i.e., cones traced by targets). In addition to the target motions, we extended the model to explain not only the horizontal offset of the vertical axis but also the tilting of the reference axis. The Gauss–Markov model with fixed constraints was adopted for the parameters to be tightly constrained, obtaining a robust and universal mathematical relationship. It minimizes the number of parameters required to estimate the IVPs, which makes the model numerically efficient and stable.

3.1. Parameterization of the IVPs

Because it is necessary to obtain the correlation information between the co-located GNSS CORS and the IVPs, we conducted a simultaneous estimation for the primary unknowns (i.e., both IVPs), coordinates of pillars, and targets based on the mathematical relationship. The ground observations were combined with two constraints: (1) the partial minimum norm least-squares solution (PMINOLESS) [25] on the pillars and (2) the fixed constraints on the target motions. The second constraint for both the horizontal and vertical targets were similarly modeled as cones (a homogeneous “conical model”).

On a preferential basis, the Gauss–Markov model with fixed constraints is represented by [22,25,26]:

$$\begin{array}{cc}\hfill \underset{n\times 1}{\mathit{y}}& \hfill =A\underset{m\times 1}{\mathit{\xi }}+\mathit{e},\mathit{e}\sim (\mathbf{0},{\sigma }_0^2{P}^{-1})\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \underset{l\times 1}{{\mathit{\kappa }}_0}& \hfill =K\underset{m\times 1}{\mathit{\xi }},\hfill \end{array}$$

(1b)

where A is the $n\times m$ design matrix, $\mathit{y}$ and $\mathit{e}$ are the observation vector, and the corresponding error vector, respectively, along with the $n\times n$ cofactor matrix $P^{-1}$ and the unitless variance component of ${\sigma }_0^2$. The $l\times m$ matrix K and the $l\times 1$ vector ${\mathit{\kappa }}_0$ are the coefficient matrix and constant vector, respectively, for the constraints from the geometric relationship of the cones. The explicit solution of the model (given in Equation (1)), along with the corresponding dispersion matrix, can be found in many references, as given by [21,25]

$$\begin{array}{cc}\hfill \widehat{\mathit{\xi }}& \hfill =N_K^{-1}\mathit{c}+N_K^{-1}{K}^T{\left(K{N}_K^{-1}{K}^T\right)}^{-1}\left({\mathit{\kappa }}_0-K{N}_K^{-1}\mathit{c}\right)\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill D\left\{\widehat{\mathit{\xi }}\right\}& \hfill ={\sigma }_0^2·\left[N_K^{-1}-N_K^{-1}{K}^T{\left(K{N}_K^{-1}{K}^T\right)}^{-1}K{N}_K^{-1}\right],\hfill \end{array}$$

(3)

where $N_K:=A^{T}PA+K^{T}K$ and $\mathit{c}:=A^{T}P\mathit{y}$. The unknown parameters $\mathit{\xi }$ can be described more clearly as follows:

$$\mathit{\xi }={\left[\begin{array}{ccc}{\mathit{\xi }}_1^T& {\mathit{\xi }}_2^T& {\mathit{\xi }}_3^T\end{array}\right]}^T,$$

(4)

where ${\mathit{\xi }}_1$ and ${\mathit{\xi }}_2$ represent the IVP parameters for the VLBI antenna and the SLR telescope, respectively, and ${\mathit{\xi }}_3$ is related to the coordinates. The VLBI parameters ${\mathit{\xi }}_1$ can be further subdivided to denote each unknown individual:

$$\begin{array}{cc}\hfill {\mathit{\xi }}_1& \hfill ={\left[\begin{array}{ccccc}{\mathit{\xi }}_{11}^T& {\mathit{\xi }}_{12}^T& {\mathit{\xi }}_{13}^T& {\mathit{\xi }}_{14}^T& {\mathit{\xi }}_{15}^T\end{array}\right]}^T\hfill \\ \hfill {\mathit{\xi }}_{11}& \hfill ={\left[\begin{array}{ccc}\alpha & \beta & \gamma \end{array}\right]}^T\hfill \\ \hfill {\mathit{\xi }}_{12}& \hfill ={\left[\begin{array}{ccc}{r}_1& \cdots & r_{n_k}\end{array}\right]}^T\hfill \\ \hfill {\mathit{\xi }}_{13}& \hfill ={\left[\begin{array}{ccc}{c}_1& \cdots & c_{n_l}\end{array}\right]}^T\hfill \\ \hfill {\mathit{\xi }}_{14}& \hfill ={\left[\begin{array}{c}h\end{array}\right]}^T\hfill \\ \hfill {\mathit{\xi }}_{15}& \hfill ={\left[\begin{array}{cccccc}{n}_x& n_y& n_z& p_1& \cdots & p_{n_k}\end{array}\right]}^T,\hfill \end{array}$$

(5)

where ${\left[\alpha \beta \gamma \right]}^T$ are the coordinates of the IVP. The elements in ${\mathit{\xi }}_{12}$ and ${\mathit{\xi }}_{13}$ (that is, $r_k$ and $c_l$) are the slant heights of the cones, which will be explained in more detail in Section 3.2. $n_k$ is the number of targets on the VLBI antenna for azimuthal (horizontal) rotation, $n_l$ is the number of targets related to the vertical motion, and h is the horizontal offset from the vertical reference axis. The unit vector ${\left[n_x{n}_y{n}_z\right]}^T$ represents the (actual) vertical axis of horizontal rotation and the (common) normal vector for the planes comprised of the rotation. $p_k$ is the distance to the horizontal plane (or base of the cone) from the apex of the cone, which is composed of the horizontal target k. The antenna-related parameters ${\mathit{\xi }}_{14}$ and ${\mathit{\xi }}_{15}$ are described, in detail, in Section 3.3. The same structure applies to the IVP parameters of the SLR telescope (${\mathit{\xi }}_2$).

The coordinate parameter ${\mathit{\xi }}_3$ denotes the increment to the initial values for the pillars (${\mathit{\xi }}_{31}$) and targets (${\mathit{\xi }}_{32}$ and ${\mathit{\xi }}_{33}$ for VLBI and SLR, respectively). As briefly explained above, the initial coordinates of the pillars and targets were calculated from the PMINOLESS approach beforehand, which were re-estimated during the IVP estimation process, with the mathematical models as constraints. The PMINOLESS approach is represented by the minimization problem [25].

$$\min \left\{{\mathit{\xi }}^{T}S\mathit{\xi }\right\}s.t.N\mathit{\xi }=\mathit{c},$$

(6)

where S is the appropriate selection matrix (a diagonal matrix with 1 for the selected elements and 0 for the rest of elements). In this experiment, the components corresponding to the pillars were selected. N and $\mathit{c}$ are the normal matrix and the corresponding vector, respectively.

The ground observations consisted of the slope distances, horizontal/vertical angles, and height differences. Therefore, the design matrix in Equation (1a) is represented as follows:

$$A=\left[\begin{array}{ccc}\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]& \left[\begin{array}{ccc}✓& ✓& 0\\ ✓& 0& ✓\\ ✓& 0& 0\end{array}\right]\end{array}\right],$$

(7)

where ✓ symbol represents the nonzero elements of the submatrices. The first row of A corresponds to ground observations related to VLBI targets, whereas the second row corresponds to the SLR targets. The third row covers the ground observations between pillars (VLBI/SLR/GNSS), including the height differences. The third column corresponds to the pillars mentioned above; thus, it is related to all types of observations. In contrast, the fourth and fifth columns of A correspond to the VLBI and SLR targets, respectively.

3.2. Fixed Constraints

As mentioned above, the three-dimensional trace of the targets on the VLBI antenna and SLR telescope can be described by two motions: azimuthal rotation (${360}^{\circ }$ a full circle) and rotation in elevation (${90}^{\circ }$ a quarter circle). These motions are similarly represented by the “sphero-conical” model [21]; however, a homogeneous conical model was applied to both motions in this study. By omitting the argument list of the coordinates ${\left[x_i{y}_i{z}_i\right]}^T$ for compactness, the conical models can be represented by two equations $R_k$ and $C_l$ [22] as below:

$$\begin{array}{cc}\hfill R_k(\alpha ,\beta ,\gamma ,r_k)& \hfill ={(x_i-\alpha )}^2+{(y_i-\beta )}^2+{(z_i-\gamma )}^2-r_k^2=0\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill C_l(\alpha ,\beta ,\gamma ,c_l)& \hfill ={(x_i-\alpha )}^2+{(y_i-\beta )}^2+{(z_i-\gamma )}^2-c_l^2=0,\hfill \end{array}$$

(8b)

where ${\left[\alpha \beta \gamma \right]}^T$ are again the coordinates of the IVP, which are the primary unknowns; and $r_k$ ($k=1\cdots n_k$) and $c_l$ ($l=1\cdots n_l$) are the slant heights of cones traced by the azimuthal rotation (target k) and the motion in elevation (target l), respectively, both from the apex at the IVP.

The whole fixed constraints can be represented by

$$K=\left[\begin{array}{ccccc}{K}_1& 0& 0& K_3& 0\\ 0& K_2& 0& 0& K_4\end{array}\right],$$

(9)

where the elements of the matrix represent the geometric relationship between the IVP and the targets for VLBI ($K_1$ and $K_3$) or SLR ($K_2$ and $K_4$). The detailed expressions for $K_1$ and $K_3$ are shown in Appendix A as Equations (A3) and (A4), respectively.

3.3. Antenna Offset and Tilting

The vertical axis of the antenna (or the telescope) may have an offset from the reference axis (i.e., a horizontal offset from the (ideal) vertical axis). Because the horizontal offset is identical for any azimuthal direction, the offset at 0${}^{\circ }$ in its own azimuth system is parameterized for both the VLBI and SLR, minimizing the number of unknowns. To estimate the horizontal offset, the mathematical model of the constraints must be modified accordingly. However, even if there is an offset on the vertical axis, it does not affect the apex ${\left[\alpha \beta \gamma \right]}^T$ of the azimuthal rotation. Therefore, the horizontal offset is exclusively related to vertical motion, and one parameter is sufficient for the entire antenna rotation.

Let ${\left[{\alpha }^{*}{\beta }^{*}\right]}^T$ be the local $(\mathit{x},\mathit{y})$ coordinates of the IVP with offset, which can be represented as follows:

$$\begin{array}{cc}\hfill {\alpha }^{*}& \hfill =\alpha +hsin\theta \hfill \\ \hfill {\beta }^{*}& \hfill =\beta +hcos\theta ,\hfill \end{array}$$

(10)

where h is the horizontal offset of the vertical axis at an azimuth angle of $0^{\circ }$, and $\theta$ represents the azimuth of the VLBI antenna (see Figure 3). Note that the bias of its azimuth system, which will be absorbed in the final Helmert transformation into the ITRF, is neglected.

$$\begin{array}{ll}{C}_l(\alpha ,\beta ,\gamma ,c_l,h)& ={(x_i-{\alpha }^{*})}^2+{(y_i-{\beta }^{*})}^2+{(z_i-\gamma )}^2-c_l^2\\ & ={(x_i-\alpha -hsin\theta )}^2+{(y_i-\beta -hcos\theta )}^2+{(z_i-\gamma )}^2-c_l^2=0,\end{array}$$

(11)

where $l=1\cdots n_l$.

However, the vertical axes of the space geodetic sensors might have tilted from the plumbline. Only the horizontal offset has been calculated by many researchers [11,27], but not surprisingly, the tilting of the axis was not modeled. Even if the vertical axis is tilted, the angle is generally very small. The horizontal targets rotate with respect to the assumed vertical axis by azimuthal rotation; thus, each horizontal plane shares the common normal vector with its plane.

The tilting of the vertical axis can be modeled using a planar equation comprising the horizontal rotation of the targets. We added three elements for the normal vector, along with the distance to the plane from the apex of the cone (i.e., IVP), to explain the tilting of the vertical axis (see Figure 4).

Assuming that a common normal vector of all planes comprised of the horizontal targets is ${\left[n_x{n}_y{n}_z\right]}^T$, then all horizontal target points can be reduced to a single plane in 3D space using the vertical distance ($p_k$) of the cone. The planar equation for each target $k(k=1\cdots n_k$) can be represented by

$$\begin{array}{cc}\hfill H_k(\alpha ,\beta ,\gamma ,n_x,n_y,n_z,p_k)& \hfill =n_x(x_i-\alpha -n_x{p}_k)\hfill \\ & \hfill +n_y(y_i-\beta -n_y{p}_k)\hfill \\ & \hfill +n_z(z_i-\gamma -n_z{p}_k)=0,\hfill \end{array}$$

(12)

where ${\left[\alpha \beta \gamma \right]}^T$ represent the coordinates of the IVP. Figure 5 shows the overall schematic diagram of vertical movement of targets at azimuth ($\theta$) under horizontal offset (h).

To secure the relationship of the horizontal planes, the normal vector should be a unit vector, as represented by

$$N(n_x,n_y,n_z)=n_x^2+n_y^2+n_z^2-1=0.$$

(13)

Since the direction of the normal vector can be either upward or downward, we applied the additional constraint of ${\left[001\right]}^T$ with appropriate stochastic uncertainty to avoid the potential singularity (i.e., the direction of the normal vector to be upward). As a result, the distance to the plane ($p_k$) is estimated to be negative when the target is located below the IVP.

The complete model for estimating the IVPs is summarized in Appendix A in linearized form.

4. Results and Discussion

4.1. Estimation of the IVP Parameters

The main concern of IVP estimation at co-located stations is obtaining the full variance-covariance matrix for the local tie vector between the geodetic sensors, along with the primary unknowns, that is, the coordinates of the IVPs. Hence, the coordinates of the pillars and targets should be estimated simultaneously with IVPs.

Table 1. Number of unknowns for each IVP, coordinates of the pillars, and targets.

Category			# Points	# Unknowns
IVP	VLBI	$\alpha ,\beta ,\gamma$	−	3
		$r_k$	−	11
		$c_l$	−	9
		h	−	1
		$n_x,n_y,n_z$	−	3
		$p_k$	−	11
		$\Sigma$	−	38
	SLR	$\alpha ,\beta ,\gamma$	−	3
		$r_k$	−	2
		$c_l$	−	6
		h	−	1
		$n_x,n_y,n_z$	−	3
		$p_k$	−	2
		$\Sigma$	−	17
Pillars		VLBI	5	15
		SLR	3	9
		GNSS	1	3
		$\Sigma$	9	27
Targets	VLBI	Horizontal	208	624
		Vertical	716	2148
		$\Sigma$	924	2772
	SLR	Horizontal	72	216
		Vertical	1201	3603
		$\Sigma$	1273	3819
Total				6673

shows the number of unknowns for each IVP, pillar, and target coordinates. A total of 55 unknowns modeled both IVPs (38 for VLBI and 17 for SLR), which are classified according to their types, that is, the slant and vertical distances of the cones and the normal vector of the horizontal planes.

Table 2. The number of required parameters with respect to the conventional 3D circle fitting type approach (k is the number of targets for AZ rotational axis, l is the number of targets for EL rotational axis, and m is the azimuth rotational sequences of EL axis).

Items	Conventional	Conical
IVP coordinates ($\alpha ,\beta ,\gamma$):
Vertical coord. ($\gamma$)	k	1
Horizontal coord. ($\alpha ,\beta$)	$2l\times m$	2
Normal vector	$3+3m$	3
Distance to plane from IVP	n/a	k
Horizontal offset	n/a	1
Radius of circle	$k+l$	n/a
Slant height of cone	n/a	$k+l$
Total	$2k+(2l+3)m+l+3$	$2k+l+7$

shows the total number of parameters required for the conventional and the newly proposed conical methods, respectively. For example, the conventional method requires up to 286 unknowns under the assumption that the same in-situ measurements presented in Table 1 are used for the VLBI IVP estimation (i.e., $k=11$, $l=9$, and $m=12$ due to the azimuth rotational sequences with 30${}^{\circ }$ increment). However, only 38 unknowns in total are necessary for the conical model (only about 13% compared to that of the conventional method). On top of the reduced parameters, the conical model is independent of the number of azimuth rotational sequences of the EL axis; we can increase the number of measurements without compromising the stability of the estimation process.

Table 3. Overall statistics of the adjustment processed in this study.

Items	Values
Observation types:
slant distance	4089
horizontal angle	3406
vertical angle	4081
height difference	8
Total unknown parameters (m)	6673
Observations (n)	11,584
Rank of A (q)	6614
Fix constraints (l)	2479
Datum deficiency ($m-\mathrm{rk}A$)	59
Model redundancy ($n-\mathrm{rk}A+l$)	7449
Convergence criterion ($|\widehat{\mathit{\xi }}|$)	5.00 mm
Convergence criterion ($|{\widehat{\mathit{\xi }}}_{\alpha ,\beta ,\gamma }|$)	0.50 mm
Number of iterations	6
Estimated variance component	5.4186
Norm of coordinate changes	0.3093

summarizes the overall statistics on the adjustment process in this study, including the number of observations and convergence criteria. The GNSS campaign was conducted on nine pillars with Trimble R8 receivers, including the GNSS (SEJN) station, on GPS week 2116 to 2117 (DOY 213 to 219) in 2020. The ground surveying (slant distances, horizontal/vertical angles, and height differences) was carried out from 28 September through 16 October using Trimble S8 total stations. All ground measurements ($N=11,584$) were used to estimate the 6673 unknowns, with convergence criteria of 5.00 mm for the norm of the correction on the entire unknowns and 0.50 mm for the primary IVPs only.

The detailed status of the iteration process is represented in

Table 4. Status of each iteration process (primary unknowns only). All units are meters except the unitless ${\widehat{\sigma }}_0^2$.

	VLBI			SLR
Iteration#	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\gamma }$	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\gamma }$	$|{\widehat{\mathit{\xi }}}_{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }}|$	${\widehat{\mathit{\sigma }}}_0^2$
Initial	37.6000	23.5000	194.0000	−55.5000	−149.0000	176.6000	0.0000	n/a
1	37.6301	23.4875	194.5947	−55.4437	−149.1018	176.4693	0.6208	5.5792
2	37.6303	23.4873	194.5953	−55.4436	−149.1019	176.4682	0.0012	5.3937
3	37.6303	23.4873	194.5953	−55.4434	−149.1019	176.4684	0.0002	5.4162
4	37.6303	23.4873	194.5953	−55.4439	−149.1019	176.4678	0.0008	5.4183
5	37.6303	23.4873	194.5953	−55.4432	−149.1018	176.4686	0.0011	5.4186
6	37.6303	23.4873	194.5953	−55.4432	−149.1020	176.4687	0.0002	5.4186

, showing the primary unknowns IVPs, norm of the parameter correction, and estimated variance component. Convergence was almost complete after two iterations, although corrections to the SLR IVP oscillated at approximately the $1$ mm level. In addition to the primary unknowns, all other parameters, including the estimated variance component, converged in a similar pattern.

Table 5. The final estimated parameters based on the Gauss–Markov model with fixed constraints and the corresponding precision. The IVP coordinates ${\left[\alpha \beta \gamma \right]}^T$ are given in the local coordinates system ($\mathit{x},\mathit{y},\mathit{z}$).

Items	VLBI (m)	SLR (m)
IVP coordinates
$\alpha$	37.6303 ± 0.0004	−55.4432 ± 0.0007
$\beta$	23.4873 ± 0.0004	−149.1020 ± 0.0008
$\gamma$	194.5953 ± 0.0007	176.4687 ± 0.0010
Slant heights of cones (H targets), $r_k$
#01	5.3402 ± 0.0006	0.6207 ± 0.0006
#02	5.3425 ± 0.0006	0.5959 ± 0.0006
#03	5.3349 ± 0.0006	n/a
#04	5.3344 ± 0.0006	n/a
#05	4.2934 ± 0.0006	n/a
#06	4.2950 ± 0.0006	n/a
#07	4.2978 ± 0.0007	n/a
#08	4.7499 ± 0.0011	n/a
#09	4.7525 ± 0.0012	n/a
#10	6.0082 ± 0.0009	n/a
#11	6.0161 ± 0.0009	n/a
Slant heights of cones (V targets), $c_l$
#01	5.3350 ± 0.0005	n/a
#02	5.3326 ± 0.0008	n/a
#03	5.3421 ± 0.0006	0.6493 ± 0.0006
#04	5.3461 ± 0.0008	0.6540 ± 0.0007
#05	4.3049 ± 0.0008	0.8032 ± 0.0006
#06	4.2936 ± 0.0004	0.7868 ± 0.0007
#07	4.2904 ± 0.0005	0.6153 ± 0.0006
#08	n/a	0.5936 ± 0.0006
#12	3.9991 ± 0.0007	n/a
#13	4.0063 ± 0.0008	n/a
Horizontal offset at azimuth $0^{\circ }$, h	0.0120 ± 0.0008	0.0033 ± 0.0005
Normal vector for horizontal planes	(no unit)	(no unit)
$n_x$	−0.0002 ± 0.0001	0.0005 ± 0.0009
$n_y$	0.0006 ± 0.0001	0.0008 ± 0.0009
$n_z$	1.0000 ± 0.0000	1.0000 ± 0.0000
Distances to planes from IVP, $p_k$
#01	1.8312 ± 0.0009	0.2065 ± 0.0008
#02	1.8287 ± 0.0009	0.2052 ± 0.0008
#03	1.8211 ± 0.0008	n/a
#04	1.8198 ± 0.0008	n/a
#05	1.8326 ± 0.0008	n/a
#06	1.8350 ± 0.0008	n/a
#07	1.8327 ± 0.0009	n/a
#08	− 0.0033 ± 0.0013	n/a
#09	− 0.0034 ± 0.0013	n/a
#10	− 3.1290 ± 0.0011	n/a
#11	− 3.1361 ± 0.0011	n/a

shows the final estimated values and corresponding precision for both the IVPs and related parameters. The slant heights were estimated with similar precision, regardless of the horizontal ($r_k$) or vertical targets ($c_l$). The horizontal offset was determined to be $12.0$ mm and $3.3$ mm for VLBI and SLR, respectively. As mentioned earlier, the horizontal offset is basically a constant value; thus, only one parameter was set up in a specific direction ($0^{\circ }$ of the azimuth angle in its own local frame) to reduce the number of parameters (or to avoid duplication).

The common normal vector of the horizontal planes for each geodetic sensor was estimated with a precision of ${10}^{-4}$ (no units). The magnitude of the unit vectors was 1.0, owing to the fixed constraints on the parameters (see Equation (13)). Parameter $p_k$ is the height of the cones that correspond to the distance to the (horizontal) plane from the apex at the IVP. As shown in the table, the distances for horizontal targets $\#08$ to $\#11$ were estimated to have negative values (marked in bold). This implies that the horizontal planes are formed in the opposite direction of the normal vector. In other words, when the antenna faces the zenith (elevation angle ${90}^{\circ }$ or zenith angle $0^{\circ }$), the corresponding targets are located under the IVP on the local vertical axis. Thus, the targets can be grouped by their heights, as seen in the VLBI target locations ($\#01-\#07$, $\#08-\#09$, and $\#10-\#11$) [22].

Finally, the simultaneously estimated coordinates of the pillars and targets are summarized in

Table 6. Estimated coordinates of IVPs, pillars, and targets with corresponding standard deviations. Only a small portion of the total targets is displayed for the sake of simplicity.

PID	x(m)	y(m)	z(m)	${\mathit{\sigma }}_{\mathit{x}}$(mm)	${\mathit{\sigma }}_{\mathit{y}}$(mm)	${\mathit{\sigma }}_{\mathit{z}}$(mm)
IVP (VLBI)	37.6303	23.4873	194.5953	0.4	0.4	0.7
IVP (SLR)	−55.4432	−149.1020	176.4687	0.7	0.8	1.0
SEJN (GNSS)	7.1613	−26.3158	181.1962	2.6	3.4	1.2
VP01	−0.0019	−0.0020	177.9101	0.4	0.5	0.4
VP02	46.5320	−0.0018	180.8347	0.4	0.5	0.4
VP03	64.9506	15.4718	180.9519	0.4	0.5	0.4
VP04	38.0174	65.5481	180.5847	0.5	0.5	0.4
VP05	1.0543	50.7583	185.5924	0.4	0.5	0.4
SP01	−65.1534	−121.5984	173.9717	0.7	0.8	0.8
SP02	−62.8026	−156.9983	179.3956	0.7	0.8	0.8
SP03	−41.1289	−142.0844	181.3180	0.6	0.8	0.8
PH01_000_90	42.0240	25.9084	196.4260	1.7	2.8	0.8
PH01_015_90	42.5024	24.6813	196.4269	1.0	2.9	0.8
PH01_030_90	42.6450	23.3734	196.4277	0.7	2.9	0.8
PH01_045_90	42.4428	22.0740	196.4284	1.0	2.7	0.8
⋮	⋮	⋮	⋮	⋮	⋮	⋮
SV08_340_10	−55.5786	−149.6105	176.7379	4.5	2.5	4.2
SV08_340_20	−55.5759	−149.6499	176.6459	4.5	1.9	4.5
SV08_340_30	−55.5741	−149.6726	176.5488	4.5	1.5	4.7

, showing only a small portion of the entire coordinates (9 pillars, 924 targets for VLBI, 1273 targets for SLR, and a total of 6618 unknowns). The precision of the estimated coordinates ranges from the half-millimeter level (pillars) to several millimeters (targets), depending on the number of observations and the geometric relationship. The initial coordinates were predetermined using PMINOLESS, with a selection matrix for the pillars (VLBI/SLR/GNSS). They were estimated again using the fixed constraints of the conical models. Notably, the GNSS CORS (SEJN) has relatively large sigma values because of the lack of ground survey measurements (and no distance observations).

4.2. Helmert Transformation to ITRF2014

The estimated IVPs and corresponding correlation information should be provided to the IERS in the global reference frame (e.g., ITRF2014). Thus, it is necessary to transform the estimated local tie vectors into ITRF2014, based on the GNSS processing results on the pillars. Considering numerical stability, it would be better to estimate the transformation parameters in the north-east-up (neu) frame by referring to the GNSS CORS as its origin.

A well-known relationship between reference frames (${\mathit{x}}_{from}={\left[\begin{array}{c}{x}_f{y}_f{z}_f\end{array}\right]}^T$ to ${\mathit{x}}_{to}={\left[\begin{array}{c}{x}_t{y}_t{z}_t\end{array}\right]}^T$) can be found in the IERS webpage and/or other textbooks in geodesy [1,22,23,28,29], whereas the velocity components are neglected in this experiment.

$$\begin{array}{c}\hfill \left[\begin{array}{c}{x}_t\\ y_t\\ z_t\end{array}\right]=\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\end{array}\right]+\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]+s\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\end{array}\right]+\left[\begin{array}{ccc}0& R_z& -R_y\\ -R_z& 0& R_x\\ R_y& -R_x& 0\end{array}\right]\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\end{array}\right],\end{array}$$

(14)

where T denotes the translation, s denotes the scale, and R denotes the rotation angle between frames.

As pointed out by [21], generally, the coordinates both from and to frames have intrinsic errors. However, the conventional least-squares adjustment only deals with errors in one frame. Thus, it is more plausible to adopt the total least-squares within the Gauss–Helmert model approach to account for the errors in both frames [22,26,30,31]. The estimated parameters and corresponding precision are listed in

Table 7. Helmert transformation parameters and summary of the estimation model (local frame to ITRF2014, both expressed in local east($\mathit{x}$)-north($\mathit{y}$)-up($\mathit{z}$) frame).

Parameter		Value	Sigma
Translation (mm)	$T_x$	1.45181	0.264
	$T_y$	−11.99716	0.286
	$T_z$	−4.22572	0.274
Rotation (arcsec)	$R_x$	−22.53129	0.782
	$R_y$	3.63510	2.245
	$R_z$	−3.83825	0.633
Scale (ppm)	s	−119.39539	3.175
Total unknown parameters (m)		7
No. observations (n)		27
Model redundancy ($\mathrm{rk}A-\mathrm{rk}B$)		20
Estimated variance component		1.19570

, and the residuals in both frames after transformation are summarized in

Table 8. Residuals in both frames after Helmert transformation (both in the local east($\mathit{x}$)-north($\mathit{y}$)-up($\mathit{z}$) frame).

		Local Frame			ITRF2014
Station		${\tilde{\mathit{e}}}_{\mathit{x}}$ (mm)	${\tilde{\mathit{e}}}_{\mathit{y}}$ (mm)	${\tilde{\mathit{e}}}_{\mathit{z}}$ (mm)	${\tilde{\mathit{e}}}_{\mathit{x}}$ (mm)	${\tilde{\mathit{e}}}_{\mathit{y}}$ (mm)	${\tilde{\mathit{e}}}_{\mathit{z}}$ (mm)
VLBI	VP01	−0.012	0.006	0.005	0.016	−0.003	−0.007
	VP02	−0.005	0.024	0.056	0.008	−0.021	−0.061
	VP03	−0.005	0.032	0.057	0.004	−0.024	−0.070
	VP04	0.006	0.021	0.027	−0.013	−0.036	−0.054
	VP05	−0.008	0.014	−0.017	0.013	−0.014	0.021
SLR	SP01	−0.390	0.341	0.677	0.118	−0.071	−0.191
	SP02	0.040	−0.908	−0.710	0.017	0.234	0.201
	SP03	0.175	1.132	0.417	−0.107	−0.301	−0.132
GNSS	SEJN	2.209	−12.118	−1.774	−0.052	0.148	0.192

.

The equation for the Helmert transformation provided by the IERS assumes very small rotation angles between the frames. However, the local azimuthal frame at Sejong station is not fully aligned to the local geodetic frame at the GNSS and is approximately rotated by $-16.5^{\circ }$ about the $\mathit{z}$-axis. Thus, the equation should not be directly applied in this case. Consequently, a preliminary rotation was used in advance. Thus, the corresponding error propagation should be applied to the variance-covariance matrix.

As a byproduct of the IVP estimation, hypothesis testing of the estimated variance component can be performed to statistically determine the existence of the horizontal offset and/or tilting of the vertical axis. Because the offset and tilting of the axis were modeled as fixed constraints, the following hypothesis can be easily applied [25,32]:

$$H_0:{\widehat{\sigma }}_1^2={\widehat{\sigma }}_2^2\mathrm{vs}.H_a:{\widehat{\sigma }}_1^2\ne {\widehat{\sigma }}_2^2,$$

(15)

where an appropriate $F_{\alpha }$ ($\alpha$: significance level) can be used as a test statistic.

5. Conclusions

A new conical model was proposed and developed in this study to estimate the IVPs of space geodetic sensors. The targets’ horizontal and vertical motions were homogeneously modeled as the bases of cones as the apex at the invariant point, resulting in a simple, robust, and immaculate model. The slant heights of the cones, except for the IVPs, were the primary parameters to be estimated within the conical model. More importantly, the horizontal offset and tilting of the vertical axis were estimated simultaneously, which makes the proposed method more consistent and applicable than the conventional approach for various IVP estimation purposes.

The estimation process began with ground surveying, followed by the IVP estimation and Helmert transformation of IVPs and the corresponding correlation information into the global reference frame. The tilting of the vertical axis was estimated by setting up a common normal vector for all planes comprising the horizontal targets. The normality of the vector was secured using a fixed constraint, and the direction of the vector was constrained stochastically to avoid singularity. Thus, regardless of the location of the targets, i.e., either above or below the IVP, the distance parameters to the planes (or bottom of the cones) were estimated with numerical stability.

The proposed model and procedure were successfully applied to the commonly used AZ-EL type antenna with in situ measurements. The test results showed that the coordinates of VLBI/SLR IVPs were estimated with the submillimeter levels of accuracy in each component. The full variance-covariance of the estimated coordinates was simultaneously determined, which makes the local tie vector combine with the ITRF solution. Furthermore, the ratio of the number of parameters in the conical model to that of the conventional one was only 13%, which is a huge advantage in terms of the numerical stability of the estimation process.

Finally, the conical model makes it possible to efficiently determine not only the IVPs of VLBI and SLR but the tilting angle of the azimuth axis and the offset between the axes with sufficient accuracy. In addition, it can provide reliable solutions by avoiding an ill-conditioned linear system caused by the limited motion of the EL axis. The applicability to other types of antennas and/or telescopes should be investigated further.