Operator-software impact in local tie networks

The operator-software impact describes the differences between results introduced by different operators using identical software packages but applying different analysis strategies to the same data. This contribution studies the operator-software impact in the framework of local tie determination, and compares two different analysis approaches. Both approaches are used in present local tie determinations and mainly differ in the consideration of the vertical deflection within the network adjustment. However, no comparison study has yet been made so far. Selecting a suitable analysis approach is interpreted as a model selection problem, which is addressed by information criteria within this investigation. A suitable model is indicated by a sufficient goodness of fit and an adequate number of model parameters. Moreover, the stiffness of the networks is evaluated by means of principal component analysis. Based on the date of a measurement campaign performed at the Geodetic Observatory Wettzell in 2021, the impact of the analysis approach on local ties is investigated. For that purpose, an innovated procedure is introduced to obtain reference points of space geodetic techniques defining the local ties. Within the procedure, the reference points are defined independently of the used reference frame, and are based on geometrical conditions. Thus, the results depend only on the estimates of the performed network adjustment and, hence, the applied network analysis approach. The comparison of the horizontal coordinates of the determined reference points shows a high agreement. The differences are less than 0.2 mm. However, the vertical components differ by more than 1 mm, and exceed the coverage of the estimated standard deviations. The main reasons for these large discrepancies are a network tilting and a network bending, which is confirmed by a residual analysis.


Introduction
In surveying engineering and large-scale metrology, network adjustment is a common analysis procedure to combine various observations, which are taken from several instrument stations but also from different instrument types. Beside the adjusted coordinates of the observed points, the network adjustment also provides the corresponding fully populated dispersion of the estimated parameters. Usually, the network adjustment only forms the basis for the subsequent specific analysis, treating the coordinates and the dispersion as incoming values. For instance, in deformation analysis the coordinates and the dispersion of several  Extended author information available on the last page of the article. epochs are evaluated and checked for congruence (Lehmann and Lösler 2017;Nowel 2020). As a further example coming from industrial metrology, in particular accelerator networks the network adjustment results are used to reconstruct and to align the path components of the accelerator beamline (Lösler et al. 2015;Manwiller 2021). At multi-technique stations such as the Geodetic Observatory Wettzell (GOW), the network adjustment results are used to obtain the reference points of the operated space geodetic techniques and to provide corresponding local tie vectors (Leinen et al. 2007;Lösler et al. 2016).
Since the network adjustment is -in most cases -independent from the specific task to be solved, in principle any software package can be used for analyzing the network. However, intensive comparison studies of different network adjustment applications have revealed discrepancies between the obtained results. Discrepancies result, for instance, from different preprocessing steps, the implemented stochastic model, or the functional model used within a specific application (Lösler and Bähr 2010;Durand et al. 2020). These software-dependent deviations are called the software effect (Durand et al. 2022). Differences in results also occur, if different analysis strategies are used within a given application. In contrast to the software effect, where the differences result from different implementations, here the selected parameters of the stochastic and the functional model within one and the same application lead to different results. This effect is known as the operatorsoftware impact (Kutterer et al. 2009), and is investigated in this contribution exemplified by the analysis of local tie surveying.
Local ties are vectors describing the relative distances and orientations between space geodetic techniques, such as the Global Navigation Satellite System (GNSS), Satellite Laser Ranging (SLR), Very Long Baseline Interferometry (VLBI), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS). The vector is known as one of the most crucial components within the combination process of space geodetic techniques realizing the global Earth-fixed geodetic reference system, and must be known with an accuracy of 1 mm to meet future requirements like a reliable monitoring of the sea level rise (Gross et al. 2009;Blewitt et al. 2010). Local ties are defined between the conventional reference points of the involved space geodetic techniques and can be derived at multitechnique stations using terrestrial measurements. The vector directly refers to a global geodetic reference frame. In classical terrestrial surveying, the instruments are leveled and relate to an astronomical frame. Wolf (1963) derives the spatial observation equations that relate the local terrestrial observations to such a global geodetic reference frame, considering the deflection of the vertical (DOV). Similar expressions can be found in numerous textbooks, such as Hofmann-Wellenhof and Moritz (2006, Ch. 5.10), Torge andMüller (2012, Ch. 6.2), andLeick et al. (2015, Ch. 4.4) to cite but a few. Kallio et al. (2022) evaluate different strategies to introduce the vertical deflections in the network adjustment. Neglecting the deflection of the vertical leads to systematic errors, when terrestrial observations obtained from, e.g., total stations are transformed into a global geodetic reference frame. For that reason, the consideration of the vertical deflection in local tie networks is recommended (Leinen et al. 2007;de Franca et al. 2022). Hereinafter, this approach is refereed to as DOV approach.
In industrial metrology, instruments need not necessarily to be leveled because the reference frame is usually an arbitrarily orientated object based reference frame. To relate arbitrarily oriented observations to an arbitrarily oriented reference frame, three coordinate components of the station position as well as three additional rotation angles must be considered. In total, each instrument position has six degrees of freedom (6DOF) (Calkins 2002). In local tie survey, this approach was applied at the Onsala Space Observatory to adjust the network observed by a non-leveled laser tracker (Lösler and Haas 2009). Later, Lösler et al. (2016) adopted the approach to directly relate terrestrial observations obtained from a total station to a global geodetic reference frame. An equivalent approach was recently applied at the Metsähovi Geodetic Research Station (Kallio et al. 2022). Hereinafter, this approach is referred to as 6DOF approach.
The mathematical model of the DOV approach and the 6DOF approach is identical. The only difference between both approaches is the handling of the vertical deflection within this mathematical model. Our main intention is to investigate the impact of the used network adjustment approach onto the resulting local ties, because both approaches are used in present local tie determinations but no comparison study has yet been made. Based on data, collected during a measurement campaign at the Geodetic Observatory Wettzell in the fall of 2021, the DOV approach and the 6DOF approach are compared to each other using the identical data set. Due to the identical mathematical model and the identical data set, differences in the results are caused by the operator-software impact, which is studied in the present contribution.
In the "Analysis procedure for local tie measurements" section, the analysis procedure is described. The model of the network adjustment is derived in the "Network adjustment" section, and the differences between the DOV and the 6DOF approach are addressed. The "Model selection" and "Principal component analysis" sections define comparison criteria to evaluate the adjustment results. Moreover, the mathematical models of the reference point determination of DORIS beacons, Synthetic Aperture Radar (SAR) corner cube reflectors, as well as telescopes used for SLR and VLBI are presented in the "Reference points" section. The proposed models are defined independently of the used reference frame, and are based on geometrical conditions. The "Case study at GOW" section deals with the data set and the analysis of the data. A brief description of the used data set is given in the "Data set" section. The data analysis is explained in the "Analysis and results" section, followed by the presentation and discussion of the results. It is shown that the results of the 6DOF approach are slightly better in comparison to the conventional DOV approach. In particular, the integration of the zenith angles is improved by the 6DOF approach, as indicated by the variance-components (VC) estimation. The comparison of the estimated horizontal coordinates of the determined reference points shows a high agreement. The differences are less than 0.2 mm. However, the vertical components differ by more than 1 mm, and exceed the coverage of the estimated standard deviations. These discrepancies result from a network titling and a network bending, and exceed the aimed accuracy goal of 1 mm in the position on a global scale. The "Summary and conclusion" section finally concludes this investigation.

Analysis procedure for local tie measurements
The determination of local tie vectors and the corresponding fully populated dispersion matrix is a complex and challenging task. As depicted in Fig. 1, the most important analysis steps are the adjustment of the terrestrial observation w.r.t. the local site network, the determination of the reference points defining the local ties between the hosted space geodetic techniques at the multi-technique station, and, finally, the transformation of the local results into a global geodetic reference frame. These final results, the estimated local tie vectors as well as the related fully populated dispersion matrix, are summarized in the Solution (Software/technique) INdependent EXchange (SINEX) format (Altamimi, 2008

Network adjustment
To combine terrestrial observation taken from different instrument stations and different instrument types in a consistent frame, a network adjustment has to be performed. As shown in Fig. 2, three differently oriented coordinate systems are to be considered within the functional model; the consistent reference frame of the network, the local system of the instrument station, as well as the local system of the target point. The datum of the network can be defined as a local geodetic xyz-frame, defined by the tangent plane of a point of tangency P 0 at the reference ellipsoid, a global Earth-fixed geodetic XY Z-reference frame, but also as an arbitrarily orientated frame. If a global reference frame is chosen, the final transformation of the local ties is omitted (Kallio et al. 2022). However, all coordinate systems can be converted into each other and, hence, do not affect the interior geometry of the network.
In this investigation, the datum of the network is defined by a local geodetic xyz-frame because the results are more descriptive and easier to be interpreted. Moreover, the local geodetic frame relates to the global geodetic frame via (Hofmann-Wellenhof and Moritz 2006, p. 209), , , Fig. 2 Relation between the conventional geodetic xyz reference frame defined by the point of tangency P 0 and the uvw coordinate frame of the terrestrial instrument. Neither the instrument nor the reflector is located at point p i or p j , respectively. Both differ by a distance h in the direction defined by g and the results can be rigorously converted. Here, ϕ 0 and λ 0 denote the ellipsoidal coordinates of the point of tangency The local coordinate system of a station is defined by the instrument axes. The u and w axes of the Cartesian system correspond to the tilting axis and the standing axis of the instrument, respectively, and are orthogonal to each other. The v-axis is orthogonal to u and w, and the origin is defined as the axes intersection point. If the instrument is leveled, the w-axis corresponds to the vertical axis defined by the local plumb line, and the instrument frame corresponds to an astronomical frame. The tilt between the instrument uv-plane and the xy-plane of the network can be parameterized by a rotation sequence defined by the angles ζ x and ζ y , and reads (Manwiller 2021) R ζ x , ζ y = In geodesy, the network points are realized by permanent ground markers or concrete pillars, and the instrument and the reflector are set up eccentrically. Therefore, the origin of the uvw-frame does not refer to the x y z T position of the consistent reference frame, and is shifted by a distance h along the vertical axis w. Using Eq.
(2), the position of the i-th instrument station expressed in its uvw-frame reads Similar to the uvw-frame of the instrument, there is a further tilted local frame having the reflector center at the target position as its origin, cf. Fig. 2. The reflector frame cancels out, if the distance h at the target position is zero. This is usually the case in industrial metrology applications (Manwiller 2021) but not in geodetic surveying (Heck et al. 1995;Jäger et al. 2017). According to Eq. (3), the j -th target position observed by the i-th instrument station is given by The coordinate differences u ij = u j − u i , v ij = v j − v i , w ij = w j − w i can readily be expressed by their corresponding polar representations. The observation equations for the slope distance, the direction, and the zenith angle are respectively. Here, the angle ζ z is a third rotation angle, and takes the non-designated frame orientations into account. The vertical axis w of an arbitrary point represented in the global geodetic XY Z-frame can be obtained from the last row in Eq. (1) (Hofmann-Wellenhof and Moritz 2006, p. 209). Inserting this axis in Eq. (1) converts the vertical axis to the local geodetic xyz-frame defined by P 0 , i.e., ⎛ Using Eq.
(2) together with Eq. (6) yields the tilting angles w.r.t. the reference ellipsoid, i.e., ζ Ell As shown by Manwiller (2021) both angles can be approximated by the central angles if the network extent is small, and the Earth is approximated by a sphere with radius R. According to Durand et al. (2022), the network extent should not exceed 1 km to obtain errors less than 0.1 mm. Equation (3) relates to the normal vector to the surface of the reference ellipsoid, which usually differs from the plumb line. Taking the deflection of the vertical into account yields where ζ DOV x , ζ DOV y are the component-wise angles of the vertical deflection.
In a least-squares adjustment ζ x , ζ y can be treated as deterministic parameters, if the vertical deflection is sufficiently known, or as parameters to be estimated. In this contribution, the first approach is refereed to as DOV approach and the second one is referred to as 6DOF approach. For the sake of completeness, it should be mentioned that ζ x , ζ y can be treated as stochastic parameters having variances σ 2 ζ x , σ 2 ζ y , respectively. In this case, both angles are introduced to the functional model as pseudo-observations as well as parameters to be estimated. However, for σ 2 ζ x,y → 0 the weighted sum of squared residuals tends to the DOV approach, and for σ 2 ζ x,y → +∞ the weighted sum of squared residuals tends to the 6DOF approach. For that reason, the DOV approach and the 6DOF approach are the two limiting extreme cases.

Model selection
A network should be parameterized by an adequate number of parameters and should have a sufficient goodness of fit. The weighted sum of squared residuals is an unqualified metric to evaluate models, because the risk of modeling the noise instead of the signal increases by an extended number of parameters, and may result in an overfitted model. Beside the goodness of fit, the model complexity must be taken into account. In statistic sciences, such a problem is known as model selection problem, and well-studied approaches are the Akaike information criterion (AIC) derived by Akaike (1974) or the Bayesian information criterion (BIC) proposed by Schwarz (1978). Both information criteria (IC) evaluate candidate models by scoring the goodness of fit against the model complexity. The preferable model among all candidate models is characterized by a simplicity of the model and by a high goodness of fit. This model is indicated by the smallest score. The basic equation reads Here, L denotes the maximized value of the likelihood function of the model, and p IC is an IC dependent complexity penalty, i.e., where p is the number of parameters to be estimated, and n is the number of observations. Equation (11a) describes a second-order improved complexity penalty that is highly recommended, if n is small w.r.t. p as shown by Burnham and Anderson (2002, p. 378). Assuming normal distributed observational errors, the log-likelihood function reads (Lehmann and Lösler 2016) where P −1 is the positive-definite dispersion matrix of the observations and v is the vector of observational residuals. Substituting Eq. (12) in Eq. (10) and excluding constant terms not affecting the scoring decision yields the fully equivalent expression where = v T Pv. AIC and BIC penalize the number of model parameters differently. By virtue of Eq. (11), BIC usually penalties complex models stronger than AIC, and tends to prefer models with less parameters. Moreover, AIC is preferred to select a model that represents the data well, while BIC is appropriate to select an explanatory (causal) model. If the true model is known and is part of the set of candidate models, BIC will select the true model with probability 1, if n → ∞. In contrast, AIC can reject the true model in favor of a better prediction as shown by Burnham and Anderson (2002, pp. 295ff); Claeskens and Hjort (2008, pp. 99ff). However, in practical applications, neither the true model is known nor the number of data increases arbitrarily, and one has to settle for a proper approximating model. In a limit range, both criteria lead to contrary decisions. For that reason, model selection methods should not be applied exclusively, and it is recommended to include additional criteria in the selection process. For a detailed analysis of AIC and BIC, the interested reader is referred to the contributions by as shown by Burnham and Anderson (2002) and Claeskens and Hjort (2008).

Principal component analysis
Beside an adequate number of model parameters and a sufficient goodness of fit, a network must be stiff. The stiffness of a network characterizes the risk that the geometry of a network deforms due to self-oscillation. A high oscillation of a network implies the presence of a socalled dominant weak-form (Jäger 1990). In the framework of deformation analysis, the stiffness is evaluated to distinguish between weak-forms and network deformations (Schmitt 1997). In optimization computations, the stiffness of a network is analyzed to refine the datum, to optimize the observation configuration, or to improve the reliability of a geodetic network (Jäger and Kaltenbach 1990;Jäger and Vogel 1990).
As shown by Jäger (1988), there is a remarkable analogy between geodetic networks and elastomechanical structures, where the eigenvalues and there corresponding eigenvectors characterize the natural frequencies and the shapes of the vibration behavior of the structure, respectively, if the system is disturbed by attacking forces. The tendency of a geodetic network to leave its geometry is, therefore, indicated by the eigensystem, when errors occur in the observations or deformations are present. Let each station be interpreted as a link of a chain, the chain becomes more flexible as the degree of freedom of each link is increased.
For that reason, treating ζ x , ζ y as parameters to be estimated reduces the stiffness of the network and increases the risk of getting a dominant weak-form.
The eigenvalues and eigenvectors are obtained by a spectral analysis of the dispersion matrix x of the estimated points x, i.e., where λ i and m i are the i-th eigenvalue and the corresponding eigenvector of x , and I denotes the identity matrix. The i-th principal component of the network reads m i √ λ i . Let the eigenvalues be sorted in descending order, so that the first principal component of the network is defined by the largest eigenvalue λ max . This principal component is called the essential component, if it dominates the full eigenspectrum of x (Jäger 1990). According to Niemeier (1982), the ratio should be less than 40% in a stiff geodetic network.

Reference points
Conventional reference points are defined from an idealized geometric structure. It is rather an exception that such defined points can be measured directly by optical or tactile measurement systems. Thus, indirect approaches are needed. In the following sections, we present general mathematical models for determining the reference points of SAR corner cube reflectors, DORIS beacons, and telescopes used for SLR and VLBI. All presented models are defined independently of the used reference frame, are based on geometrical conditions, and are refereed to conventional reference points.

SAR corner cube reflector
Trihedral corner cube reflectors are used to evaluate the pixel localization accuracy of a SAR system. The corner cube consists of three orthogonal sides and the reference point p SAR is defined as the intersection point of these sides, which is equal to the position of the electric phase center (Balss et al. 2018). In order to validate the SAR accuracy, the position of a corner cube obtained by SAR is compared to the reference point determined by precise terrestrial geodetic surveying (Gisinger et al. 2015). Figure 3 depicts a schematic representation of a trihedral corner cube reflector. Each side of the corner cube is parametrized by a plane. The j -th plane is defined by a normal vector n and a form parameter d, i.e., Here, p is an arbitrary point lying in the plane. Depending on the used instrument type, the surface of a plane is measured either tactile using a reflector, or target-less by direct measurements. If the plane is observed tactilely by, e.g., a laser tracker, the known physical distance r Ref between the center of the reflector and the plane has to be taken into account, cf. Fig. 3. For target-less measurement systems such as laser scanners or photogrammetric systems r Ref = 0 is obsolete. The parameter d and the three components of the normal vector n are linear dependent. A common condition to prevent the ambiguity between the parameters is to normalize the normal vector where · 2 symbolizes the Euclidean norm. In this case, d represents the shortest distance between the plane and the origin, and Eq. (16) is known as the Hesse normal form. The reference point p SAR is the inner corner and, hence, is the only point lying in each plane. For that reason, the point must fulfill the condition for j = {I, II, III}. Equations (16), (17), and (18) describe the model of a three-sided pyramid, where the apex corresponds to the reference point. However, this derivation does not ensure an orthonormal basis of the normal vectors of the planes. To obtain a right-angled three-sided pyramid, three further conditions must be introduced, i.e., with i = {I, II}, k = {II, III}, and i < k.
The orientation vector o SAR of the corner cube reflector is readily obtained from the three normal vectors n j of the planes and reads In order to obtain the 12 parameters of the three planes as well as the reference point position of the SAR reflector, a generalized least-squares model with additional unknown parameters in the condition equations is required.

DORIS beacon
The conventional reference point is a non-materialized point inside the structure of the DORIS beacon. The beacon has a cylindrical form, and the reference point p DORIS is defined as the point lying on the main axis of the cylinder, which is located h c = 390 mm above a reference plane . Figure 4 depicts a schematic representation of a DORIS beacon. A red ring at the beacon surface marks a position close to the reference point, but the center of the ring does not necessarily coincide with the reference point . The phase center p 2GHz of the beacon is 877 mm above the reference plane, and is located on the main cylinder axis. Scheme of a DORIS beacon. The phase center p 2GHz is symbolized by a green star. The reference point p DORIS is the location point of the cylinder, and is depicted by a red dot. The distance between the p DORIS and the reference plane is h c = 390 mm. The orange normal vector of the reference plane is identical to the main axis of the cylinder. The radius of the cylinder is r c . At the reference plane, a terrestrial reflector with radius r Ref is symbolized DORIS beacons provide a mounting support for conventional terrestrial reflectors below the reference plane. This position is centered on the cylinder axis. If the pointing direction of the axis is known or if it is strictly vertical, this position provides an easy access to the reference point ). However, the position and the tilt of the beacon can change over time. For instance, a tilt of 0.5 • shifts the phase center by about 7 mm. For that reason, a unified approach is suggested that provides both, the reference point and the pointing direction of the axis. This approach can be easily adapted to the reference point determination of GNSS antennas.
According to Lösler (2020) the intrinsic parameters of a circular cylinder are the normal vector n c , and the radius r c . Let p c be the location position, the common equation of a circular cylinder reads where p denotes an arbitrary surface point. Due to the ambiguity of the normal vector, the vector is normalized using Eq. (17), i.e., n c 2 = 1.
The reference plane is orthogonal to the cylinder axis, so that the normal vector is identical. According to Eq. (16), the plane reads where d c is the shortest distance between the plane and the origin, and r Ref compensates for a known reflector offset, if tactile measurement systems are used. The location position p c can be chosen arbitrarily. Considering the distance h c between the plane and the reference point yields a distinct location position, which is identical to the reference point, i.e., where p c = p DORIS and h c = 390 mm. The parameters to be estimated are the position of the reference point, the normal vector parameterizing the cylinder axis, the radius of the cylinder, and the shortest distance between the plane and the origin. To estimate these eight parameters, a generalized adjustment model with conditions is recommended.

SLR/VLBI reference point
The reference point of swivel mounted telescopes used for SLR or VLBI is defined as the orthogonal projection of the secondary axis onto the primary axis of the telescope. Without loss of generality, let the primary axis and the secondary axis be the azimuth axis and the elevation axis, respectively. According to this geometrical definition, the reference point is independent from the pointing direction of the telescope. Since the reference point is enclosed by the telescope structure, it is generally not materialized and cannot be measured directly by means of, e.g., optical or tactile measurement systems (Leinen et al. 2007;Gong et al. 2014). For that reason, the reference point is indirectly obtained by analyzing the trajectories of mounted markers at the turnable part of the telescope structure, which result from the rotation of the telescope in varying azimuth and elevation orientations.
A common approach to determine the reference point is derived by Lösler (2009). Here, the reference point p Tel results from a transformation of the telescope-fixed frame to the consistent Earth-fixed reference frame, cf. Fig. 5. The mathematical model reads Here, the angles α and β compensate for the tilt between the azimuth axis and the z-axis of the reference frame, and γ is the non-orthogonality between both telescope axes.
The k-th telescope pointing direction is defined by the azimuth angle κ k and the elevation angle ω k . The angle κ 0 describes the unknown orientation between both frames, and κ * k = κ k + κ 0 is the oriented azimuth angle w.r.t. the Earth-fixed frame. Rotation sequences are parametrized by basic rotation matrices R. The axis-offset is e AO . Vector q i Fig. 5 Geometric relations between the telescope-fixed frame and the Earth-fixed reference frame, which is used to determine the reference point p Tel of SLR/VLBI telescopes. The swivel-mounted telescope rotates around the azimuth and elevation axis with angles κ and ω, respectively, and forms the trajectory of a mounted marker p i . The axis-offset is e AO is the position of the i-th marker defined in the telescopefixed frame, and p i,k is the corresponding position in the Earth-fixed frame. The coordinates of the reference point as well as the nuisance parameters can be adjusted by treating the marker positions as observations. Initially, the model was derived for reference point determination at VLBI radio telescopes (Lösler 2008). However, the model is also valid for telescopes used for SLR as recently demonstrated by Lösler et al. (2018Lösler et al. ( , 2021 and Eschelbach and Lösler (2022). In contrast to geometrical approaches like the circle-fitting approach (Leinen et al. 2007;Gong et al. 2014), Eq. (25) does not require predefined telescope orientations, and is suitable for an automated reference point determination during the regular station process (Ning et al. 2015;Lösler et al. 2016). Moreover, the model is not restricted to a specific telescope type and contains important additional parameters like the axis-offset e AO , which, if unconsidered, biases the global position (Combrinck and Merry 1997), or the zero-pointing direction κ 0 , which is used to reduce the pointing error of the telescope (Zhang et al. 2021).

Sequential estimation of the parameters
Local ties are introduced to the combination process to overcome the weak physical connection between different space geodetic techniques. For a rigorous combination, the dispersion of the local tie vectors must be considered (Altamimi et al. 2016). In the previous section, we provide a general framework for deriving conventional reference points of different space geodetic techniques. The estimates of the network adjustment, x and x , are treated as incoming values. It is assumed that x and x are reduced to quantities needed for reference point determinations. In fact, determining the reference points individually yields the dispersion of the reference points but mutual correlations between the reference points get lost, and the dispersion of the resulting local tie vectors becomes a block-diagonal structure. By determining all reference points in a single adjustment, the correlations implied by the dispersion matrix are rigorously propagated (Ma et al. 2018). However, such a model strongly depends on the number of hosted space geodetic techniques, and has to take site specific strangeness into account. Considering all contingencies is challenging or nearly impossible. For that reason, a more flexible and adaptable procedure is suggested.
The general adjustment model with additional parameters in the condition equations results from the method of Lagrange multipliers given by Here, denotes the weighted sum of squared residuals, x is the vector of observations, and β T = β T 1 β T 2 are the parameters to be estimated. The target function is combined with the implicit functional model f and the condition equation c using the Lagrange multipliers λ T = λ T 1 λ T 2 . Whereas β 1 appears in f and c, β 2 denotes the additional parameters, which appear only in c. To obtain the minimum, the partial derivations w.r.t. β, x, and λ are evaluated and equated to zero. The resulting normal equation system reads (Mikhail 1976, p. 234 Here, matrices A and B consist of the partial derivation of f w.r.t. to the unknown parameters β 1 and the observations x, respectively, and matrix C consists of the partial derivation of c w.r.t. the parameters β. The vector of increments is denoted by β . Vectors w 1 = f − Bv x and w 2 contain the misclosures w.r.t. f and c, respectively. Equation (27) is solved iteratively by stripping the new iterates of their randomness, and treating them as improved approximations into the next iteration step, i.e., wherev x = x B Tλ 1 . The final step yieldsβ,x, andλ, as well as the related dispersion matrices, β , x , and λ . It is strongly recommended to evaluate the strength of the introduced constraints using, for instance, hypothesis tests as derived by Lehmann and Neitzel (2013).
In order to obtain the stochastic dependencies betweenβ and x, Eq. (27) is rearranged and is expressed as a linear equation system that only depends on the observations x (Mikhail 1976, p. 116f), i.e., where the vector m contains the deterministic parts of the model, and M = B x B T . Matrices E are selection matrices and remove quantities, which become meaningless for further analysis steps like the nuisance parameters. The vectorβ = Eββ contains the coordinates of the reference point, and the vectorx = Exx denotes the remaining points of the network adjustment. Applying the propagation of uncertainty yields where the block-matrices are respectively.
According to Eqs. (29), (30) a certain reference point is estimated and rigorously combined with the remaining points of the network taking mutual correlations into account. This procedure is repeated until all reference points are estimated by setting x ← y and x ← y . At the end of this procedure, vector y represents the local ties and y is the corresponding fully populated dispersion matrix. This procedure is independent of the computation order of the reference points.

Case study at GOW
In order to evaluate the impact of the used network adjustment approach onto the reference points, the data set of the 2021 local tie campaign carried out at the Geodetic Observatory Wettzell (GOW) is chosen. GOW is a multitechnique station and operates instruments of all basic space geodetic techniques (Hugentobler et al. 2011). The local site network depicted in Fig. 6 consists of about 30 concrete pillars, surrounding the hosted space geodetic techniques, that are one DORIS beacon (WEUC), two SLR telescopes (SOS-W, WLRS), three VLBI radio telescopes (RTW, TTW-1, TTW-2), and several GNSS antennas (WTZA, WTZR, WTZS, WTZZ, WLRG). Moreover, GOW supports SAR missions by two corner cube reflectors (SAR1, SAR2). In this investigation, we restrict ourselves to the analysis of the DORIS beacon WEUC, the corner cube reflector SAR2, the radio telescope TTW-2, as well as the SLR telescopes SOS-W and WLRS, because indirect measurements were performed for these techniques, which allow the determination of conventional reference points by means of datum-independent geometrical conditions.

Data set
The measurement campaign was carried out in the fall of 2021. The survey data form the basis for the local tie vectors for the upcoming International Terrestrial Reference Frame. The measurements were carried out using the Leica total stations TS50 and TS60. Both instruments are specified identically by 0.6 mm + 1 mm km −1 and 0.5 for distance and angle measurements, respectively (Leica 2013;2015). Moreover, the novel multi-lateration system DistriMetre was used for data acquiring at GOW for the first time. This instrument is a prototype developed by CNAM. In controlled environments, the uncertainty of the distance measurement is about 5 um (Guillory et al. 2020b;. However, under the environmental conditions of the measurement campaign at GOW, the uncertainties on the measured distances were larger. They were mainly affected by the limited knowledge of the air temperature used for the air refractive index determination. The multi-lateration was applied to improve the reference point determination by high-precise distance measurements at the southern radio telescope TTW-2. The remaining reference points as well as the pillar network were observed by the total stations. All pillars were used as instrument station except for those equipped with GNSS antennas. As already mentioned, the estimation of the extra rotation parameters requires an observation configuration that is sensitive for these parameters within the 6DOF approach. For that reason, measurements were taken from each pillar to at least three other pillars. To ensure a reliable integration of each instrument station, measurements were carried out to all visible network pillars. Figure 7 depicts the realized observation configuration of the GOW pillar network. In total, more than 31,000 terrestrial observations were conducted. The International Association of Geodesy (IAG) recommends the dispersion model derived by Ciddor (1996Ciddor ( , 2002 and Ciddor and Hill (1999) for high-precise geodetic applications to compute the first velocity correction (IAG 1999). Following this recommendation, temperature, pressure and humidity were measured at each instrument station using a GFTB 100 (Greisinger), because these meteorological parameters are the most crucial parameters affecting the first velocity correction (Rüeger 1996, Ch. 5.8). Due to its smaller impact on the first velocity correction, carbon dioxide was measured stationary close to pillar 34 using a CO 2 -Meter GC-2028 (Lutron). Table 1 summarizes the manufacturer specifications of the meteorological measurement sensors. Using these meteorological parameters, the distances of the total stations were corrected using the procedure given by Pollinger (2020). Furthermore, all measurements of the total stations were carried out in both faces, to reduce systematical deviations caused by mechanical deficiencies of the total stations. For data analysis mean values from repeated measurements of full sets of polar observations were calculated and introduced to the network adjustment.
The distances of the DistriMetre were corrected for mechanical misalignment and systematical errors using apriori values taken from laboratory calibrations. The group refractive index of air was calculated with the modified Edlén's formulae (Bönsch and Potulski 1998), using the recorded environmental parameters. A detailed description of the multi-lateration measurement system DistriMetre and an evaluation of the uncertainties are given in the contributions by Guillory et al. (2020a, b).
In total, the preanalyzed data set consists of 1310 points which were connected by n = 7383 corrected and averaged terrestrial observations. These observations were introduced to the network adjustment.

Analysis and results
The network adjustment is realized by the in-house software package Java·Applied·Geodesy·3D (JAG3D 2022) using the functional model described in the "Analysis procedure for local tie measurements" section. The datum of the network adjustment is defined by the tangent plane of P 0 , which is equal to the centroid of the closed polygon defined by the outer vertices of the network. In Fig. 6, a green hexagon symbolizes the position of the principal point P 0 . The concrete pillars of the network realizing this datum within the free network adjustment are depicted by red triangles. The reference ellipsoid is GRS80. According to Kallio et al. (2022), the deflections of the vertical are ζ DOV x = −5.85 and ζ DOV y = +2.18 , and a uniform tilt of the points of the entire network is assumed. It must be mentioned that uniformly applied vertical deflections do not effect the interior network geometry but only the orientation of the network. By applying the DOV approach, these vertical deflections are assumed to be deterministic. In the 6DOF approach, the vertical deflections are additional parameters to be estimated.

Network adjustment
In order to study the impact of the functional model on the estimates, the datum of the network, the observations, and  Table 2 Components of the stochastic model used for distances s, directions τ , and zenith angles ϑ observed by the total stations (TS) and the DistriMetre (DM) system, respectively. The zero point offset uncertainty is denoted by σ a , and σ b is the distance dependent component where s is the slope distance, and ρ converts the angle unit from radians to arc seconds. Numerical values for σ a and σ b are summarized in Table 2. During the network adjustment, the components of the stochastic model are evaluated by means of variancecomponents (VC) estimation (Crocetto et al. 2000). Table 3 Result of the VC estimation. For each approach, the number of observations n, the effective number of observations n eff , the redundancy r, the weighted sum of squared residuals , and the resulting variance-componentσ 2 = /r are given component-wise because small deviations do not significantly affect the estimates and are tolerable (Durand et al. 2022). Table 3 summarizes the total number of observations n, the total redundancy r, the weighted sum of squared residuals , and the estimated VC. Points used to determine reference points are generally not redundantly observed and therefore do not contribute to the weighted sum of squared residuals . By excluding these observations, remains unchanged, and the effective number n eff of observations is obtained. Table 4 gives an instrument-related overview. The estimated VC satisfy Eq. (33) and, with the exception of the VC of the zenith angles, are below the expected value E σ 2 = 1.0, which means that the stochastic model is slightly too pessimistic. However, it should be noted that the VC estimation is based only on the dispersion of the observational residuals, which are not sensitive to all measurement discrepancies, as recently discussed by Neitzel et al. (2022). For that reason, a further adaption is omitted. By comparing the component-wise estimation between both approaches, the differences in the distance VC and direction VC are negligible. Noticeable differences can be found for the VC of the zenith angles, which are most affected by the handling of the vertical deflections.
The mean standard deviations of both approaches are summarized in Table 5. Whereas the horizontal components yield almost comparable results, the vertical component is larger by about 0.1 mm. Moreover, the DOV approach leads to slightly smaller standard deviations than the 6DOF approach, but both approaches fulfill the accuracy requirement of <1 mm. Figure 8b depicts the estimated vertical deflections obtained from the 6DOF approach. For most of the vertical deflection vectors, the orientation and the length match very well in comparison to the gravimetrically obtained values shown in Fig. 8a. The reliability of the estimated vertical deflections depends on the measurement configuration. Pillars, which are connected to a large number of pillars, yield a better match than less observed pillars located, for instance, at the boundary of the network as indicated by orange rings in Fig. 8b. It should be noted that the estimated vertical deflections are a combination of overlaying effects. The standing axis error, for instance, describes an additional deviation from the local plumb line, and is not compensated by measuring in both faces. Moreover, the orientation and the length of the estimated vertical deflections depend on the introduced approximation coordinates and, therefore, on the network datum. From a mathematical point of view, the DOV approach neglects these additional systematical deviations and, thus, is less rigorous.
The estimated vertical deflection vectors show a slight torsion, which may be composed of a network tilting and network bending. A network tilting occurs, if an inappropriate datum is specified. A network bending, on the other hand, results from an improper functional model that, for instance, forces the interior geometry of the network. A spatial six parameter transformation is performed to separate both effects, considering the estimated datum points of the 6DOF and the DOV approach as source and target system, respectively. The estimated tilting angles are ζ x = 1.27 and ζ y = −0.67 , and result in a vertical discrepancy of about 2 mm at a distance of 300 m. The remaining residuals of the points confirm the assumption of an additional network bending. The residuals of the horizontal components are negligible. The residuals of the vertical component vary in a range of about ±0.8 mm and are depicted as contour plot in Fig. 9. Smallest values can be found in the center of the network close to P 0 . The residuals increase with increasing distance from P 0 . Such a bending may be caused by a dominant weak-form indicated by the first principal component.

Model selection and stiffness
On the one hand, the 6DOF approach reduces the weighed sum of squared residuals by treating the vertical deflections as additional parameters to be estimated. On the other hand, using a more complex model increases the risk to get an overparameterized model. Based on AIC and BIC, the goodness of fit is evaluated w.r.t. to the model complexity using Eq. (13). Moreover, the stiffness of the network is affected, if the number of model parameters is changed. According to Eq. (15), the stiffness is characterized by the largest eigenvalue w.r.t. the full eigenspectrum. Table 6 summarizes the score values of the model selection methods as well as the analysis result of the first principal component. AIC, BIC and pc max are given for the complete network. Here, AIC and BIC select the 6DOF approach, indicated by smaller IC values. By virtue of Eq. (13), the evaluation of the goodness of fit depends on and n as indicated by the first term. Non-redundant observations do not contribute to but decrease the first term and, thus, the score value. For that reason, Table 6 contains also the score values calculated from the redundantly observed network. Again, AIC favors the 6DOF approach, but BIC selects the DOV approach. The AIC is recommended for prediction and, thus, the observations are better represented by the 6DOF approach. However, this approach is not necessarily closest to the (unknowable) true model as assigned by the BIC. The analysis of the first principal component pc max shows a reduction in the stiffness, if the vertical deflections are estimated within the network adjustment. This becomes even clearer, when the redundantly observed network is evaluated. The ratios pc max increase by a factor of 1.4 and 2.4, respectively, but remain well below 40%. Neither approach has an essential first principal component, and the network bending depicted in Fig. 9 does not result from a dominant weak-form. Since only two approaches are compared, the cause cannot be clearly determined, and   . 9 Contour plot of the vertical residuals obtained from a spatial transformation of the corresponding datum points of both approaches needs further investigations. However, these differences result from the operator-software impact and occur, if different operators apply different analysis strategies to the same data using identical software packages.

Reference point determination
Based on the results of the network adjustment, the reference points are determined using the suggested sequential analysis procedure. The estimates obtained from the DOV approach are summarized in Table 7 w.r.t. the coordinates of SAR2. Additionally, the coordinate differences, = 6DOF − DOV, and the standard deviations resulting from the 6DOF approach are given. The differences of the horizontal coordinate components are quite small and do not exceed 0.2 mm. However, the vertical coordinate components differ significantly. By far the largest discrepancy of about 2.5 mm can be found for the DORIS beacon WEUC, and exceeds the coverage of the estimated standard deviation of 0.4 mm. This discrepancy results from the pillars 25 and 26 in the immediate vicinity to WEUC, which show similar discrepancies in the vertical components and have been used to observe the surface of the beacon. For the SLR and VLBI telescopes as well as the corner cube reflector SAR2, the horizontal components are almost identical and the vertical components differ by about ±0.5 mm. Both, pillars and obtained reference points, are affected by the network bending (see Fig. 9) and the network tilting.
The deformation of the interior geometry of the local ties is evaluated by the residuals of a spatial six parameter transformation. Table 8 contains the component-wise residuals. Largest values can be found for the vertical component. All residuals are less than 0.5 mm. Thus, the bending is small w.r.t. the total discrepancies given in Table 7, and the differences between both approaches mainly result from the network tilting. The DOV approach transforms the network by four parameters onto the approximation values defining the datum of the network, namely three translation parameters as well as one rotation parameter. In contrast, the 6DOF approach considers two more rotation parameters, and introduces six parameters. Due to the extra parameters, the network fits better to the defined datum, if reliable observations are collected that are sensitive for the extra rotation parameters. However, if the datum definition is weak, the network orientation is also weak. Due to terrestrial measurements being insensitive to define an absolute frame, the datum itself can only be evaluated in a global context.
The fully populated dispersion matrix of the local ties results from the suggested sequential analysis procedure. This matrix describes the uncertainties and the linear dependencies of the estimated parameters. Both considered approaches yield comparable correlations between the estimated reference points, as depicted in Fig. 10. About 96% of the correlations are less than 15%. Larger dependencies arise between the coordinate components of one and the same reference point but largely decay between different reference points. The maximum correlation is 33% and 43% for the DOV and the 6DOF approach, respectively. These minor dependencies result from the measurement configuration at GOW. The reference points are observed independently from each other using the concrete pillars in the immediate vicinity to the particular space geodetic technique. Simultaneous measurements to several space techniques are generally not possible, resulting in only minor dependencies.
Even though both approaches yield reliable results, the numerical differences are noticeable and can be addressed to the operator-software impact. The variations in the results of different analysis strategies can by far exceed the measurement uncertainties specified in the stochastic model during the adjustment process. In the framework of local tie determination, such discrepancies can be larger than the aimed accuracy requirements and induce errors in the combination process.

Summary and conclusion
Most applications in engineering geodesy and industrial metrology perform network adjustments. A network adjustment combines the collected observations w.r.t. the relating measurement uncertainties and yields a consistent Cartesian frame. As shown by Durand et al. (2020), the results can differ, if the same data are analyzed using different adjustment software packages. These deviations are known as software effect and result from different implemented functional or stochastic models (Durand et al. 2022). Even the use of a certain software package leads to different results, if different operators apply different analysis strategies to the same data. This effect is known as operator-software impact (Kutterer et al. 2009). In this investigation, the operator-software impact was investigated in the framework of local tie determination for the first time. For this purpose, a universal procedure to obtain reference points was presented. All suggested models within the procedure are defined independently of the used reference frame, and are based on geometric conditions. The estimates of the network adjustment are treated as incoming values. The procedure is sequentially applicable and yields the local ties as well as the fully populated dispersion matrix. Two different analysis strategies were compared, which mainly differ in the consideration of the vertical deflection within the network adjustment. The mathematical model of both approaches is identical. Whereas the DOV approach treats the deflections of the vertical as deterministic parameters, the 6DOF approach considers the vertical deflections as additional parameters to be estimated. Increasing the number of model parameters within a more complex model increases the risk to get an overparameterized model as well as to obtain a dominant weak-form. Both risks were addressed using information criteria and principal component analysis.
To evaluate the impact of the analysis strategy on the reference points of space geodetic techniques, both approaches are applied to the identical data set of a local tie campaign performed at the Geodetic Observatory Wettzell in 2021. The network was analyzed using the in-house software package JAG3D. The information criteria AIC and BIC evaluate both approaches differently. However, one can conclude that the observations are better represented by the 6DOF approach, but this approach is not closest to the true model. The estimated vertical deflections of the 6DOF approach agree well with gravimetrically obtained values. Moreover, neither approach has an essential first principal component, and both approaches yield stiff networks.  The comparison of both approaches reveals discrepancies in the determined reference point coordinates. Whereas the horizontal coordinate components are almost identical, the vertical coordinate components deviate by more than 1 mm, and exceed the coverage of the estimated standard deviations. These discrepancies result from a network tilting and a network bending. A network tilting occurs, if an inappropriate datum is specified. Due to the extra parameters of the 6DOF approach, the network fits better to the defined datum. However, if the datum definition is weak, the network orientation is also weak, and network tilting occurs. Unfortunately, the datum can only be evaluated in a global context because terrestrial measurements are insensitive to define an absolute frame. A network bending, on the other hand, results from an improper functional model that forces the interior geometry of the network. It is well-known that deterministic model parameters affect the interior network geometry. However, treating these parameters as additional unknowns requires reliable observations that are sensitive to the parameters. The network bending of the network at the Geodetic Observatory Wettzell is less than 0.5 mm and, thus, both approaches yield almost comparable results and fulfill the aimed accuracy goal. The crucial part of the discrepancies results from the datum definition.
A general recommendation is challenging, because the measurement data allows a limited number of meaningful approaches. Further investigations and comparisons are needed to provide suitable analysis procedures and guidelines, which reduce the datum problem in local tie networks, and the discrepancies caused by the operator-software impact. Future work will focus on the benefit of introducing additional leveling data to the adjustment process, which could improve the estimated extra rotational parameters and, therefore, the stiffness of the network. Including leveling data could lead to a clear information criterion decision and to reduce the discrepancies between both approaches under investigation. Moreover, the dependency between the adjustment approach and the network extent as well as the network configuration should be addressed.
use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.