A Method for the Precise Coordinate Determination of an Inaccessible Location

Surveyors are occasionally tasked to with determining the coordinates of inaccessible locations or points in civil engineering applications, ground control points for photogrammetry or LiDAR data acquisition, among others. The present work outlines and investigates a novel method for estimating the GNSS coordinates of an inaccessible location where a surveying instrument cannot be set up. The procedure is based on the well-known surveying intersection method and data extracted from an Earth Gravity Model (e.g., EGM 2008). The location’s coordinates are obtained from the least-squares adjustment of the angles and distances measured from at least two sites to the unknown point using a total station, within the framework of the Gauss–Helmert method. Field tests confirmed that the accuracy of the determined coordinates of the inaccessible point is at the level of 1 cm. The proposed method bypasses standard coordinate transformation steps performed with the traditional approach, directly producing geocentric coordinates of the unknown points.


Introduction
Determination of the three-dimensional coordinates of locations or unknown points is a fundamental task in professional surveying. Various measurement techniques can be used for this purpose, among which the GNSS (Global Navigation Satellite System) method is probably the most obvious. However, the GNSS method cannot be used if the GNSS receiver is obstructed (e.g., by tall buildings in the city canyons) or if the location is inaccessible. Nevertheless, the goal is still achievable if the inaccessible point is visible from at least two points using a GNSS station, and a total station is on hand to measure angles and distances to the unknown point from the GNSS accessible points.
The problem of integrating various geodetic data types has been an active research topic for a long time. Some relevant publications on this topic include, e.g., [1][2][3][4][5]. These publications demonstrate different approaches to combining various spatial data types with joint adjustment and visualisation. The authors analysed different GNSS networks, which were a base for other measurements (for example, total station) and for capturing inaccessible points. The results are also used in engineering tasks (like deformation monitoring) or for providing control surveys (e.g., rock block displacements) in the mountains. Data integration and fusion procedures were presented regarding integrating available geoid models. Then, the objects were spatially analysed and adjusted in a coherent, unique reference system, giving more satisfactory results than by processing the data separately.
Implementing combined surveying methods, especially total station and GNSS methods requires continuous scientific research and validation under field conditions. Moreover, similar approaches can successfully be utilised in the mining industry [6]. Examples of compelling applications and reports on developing spatial data integration methodologies are published in [7]-both cases concern deformation surveys and precise levelling based on the joint adjustment of multi-source data. In the case of open pit mines, one can observe a significant influence of atmospheric refraction and, at long distances, the Earth's curvature (plumb-line deviations). Similar problems are also visible when performing levelling across wide rivers. In this case, geodetic works must establish special control networks on two opposite riverbanks and utilise unique instrumental designs (e.g., integrated prisms aided by LED lighting). Today, such measurements can be performed with millimetre-level accuracies.
However, such cases require a skilful fusion of the mentioned aspects, supported by advanced computational methods. In addition to GNSS and total station measurements, other technologies like laser scanning or photogrammetry are also worth noting. Some examples of utilising the above surveys are presented in [8,9]. The first paper focuses on the issue of inventory surveying of façade scaffolding. The authors considered various land surveying methods such as photogrammetry, laser scanning or total station. The developed method for integrating the abovementioned data allows for the effective monitoring of scaffolding on construction sites. In the second article, the authors demonstrate their newly developed calculation and data adjustment method for precise object dimensioning by employing close-range photogrammetry. This solution utilises a modified Gauss-Markov approach to adjust the spatial control network deployed directly on the surveyed object. Both demonstrated cases offer singular millimetre accuracies. In modern engineering geodesy, the development of computational methods-especially planimetric and numerical-has become particularly important.
The spatial intersection method, well known to surveyors, allows for the determination of three-dimensional point coordinates based on angle and distance measurements from at least two points where total station/GNSS receivers can be placed (e.g., [10][11][12][13]). In addition, in [14], the authors examined the accuracy of coordinates of inaccessible points determined using the intersection method. Such an approach can be successfully employed where measuring direct distances is impossible. Hence, it finds numerous applications in civil engineering tasks, including the setting out of control networks, control surveys or similar jobs [15]. In [16], however, a novel concept in modern geodesy, the Gauss-Helmert model (GHM), was used. The GHM [17] is a variation of the least-squares method designed to minimise errors in all variables, including parameters and observations. It is formulated as linear conditional equations with unknown parameters. The proposed intersection method for the coordinate determination of an inaccessible location uses linear conditional equations. Hence, the GHM is a relevant approach to solving the system of observation equations.
Textbooks provide in-depth background information on this approach, e.g., [18]. In addition, the properties of this method have also been analysed in the literature, e.g., [19][20][21][22]. Various authors have also studied related issues, including similarity transformation problems under the GHM, proposed alternative solutions [23,24], and compared with other similar approaches, e.g., Gauss-Newton [25]. Also, some more advanced topics, including the GHM with a singular dispersion matrix [26], parameter estimation within the GHM using the Pareto optimality [27] or robust parameter estimation of the non-linear GHM [28], have been investigated.
The present paper introduces the GHM for the intersection method. The proposed approach integrates GNSS observations and the total station's deflections of the vertical axis, determined based on a global Earth Gravity Model (EGM) [29]. The authors' primary motivation is to derive the most general GHM of observational equations binding various possible measurements in the proposed spatial intersection. Moreover, one of the leading practical aspects was the ability to offer an efficient method for surveying inaccessible points in measuring services employing real-time GNSS receivers. A similar approach was already tested and published in [30]; however, the demonstrated technique has been completed using the GHM. Above all, it should be noted that the classic surveys, together with their modern forms using, e.g., aerial or terrestrial laser scanning, can be reliably integrated and combined with GNSS by expressing the results in a common geocentric datum. Thus, we introduced the geocentric intersection method, which includes rigorous accuracy assessment.
The specific processing steps that are performed in the present paper can be synthesized as follows:

•
The problem: to obtain precise geocentric coordinates of inaccessible locations for photogrammetry or LiDAR survey. • A method is proposed based on the directional intersection of the inaccessible locations using a total station from two stations with measured GNSS coordinates.

Basic Relationships
Using a total station, one can measure the horizontal α and vertical β angles and distance s from point P to target Q ( Figure 1). Assuming that the vertical angle β is corrected due to the atmospheric refraction and the instrument/mirror height (i)/(j) is considered, the coordinates (x, y, z) of target Q are computed from (1): The coordinates x of target Q are referenced to a local coordinate system (x, y, z). They must be transformed into the destination reference system. As a destination reference system, we use the geocentric World Geodetic System 1984 (WGS84), which the GNSS uses.
Referring to Figure 1, the origin of the local reference system is located at P. Point P with geocentric coordinates (X, Y, Z) in the GNSS reference system and the geodetic latitude ϕ and longitude λ. The orientation of the local system (x, y, z) in the geocentric system (X, Y, Z) is defined by the three orientation angles: Σ, ξ, η. The so-called total station horizontal orientation constant is represented by the angle Σ. The ξ and η are vertical deflection angle components (north/east) from the normal to the ellipsoid at P.
The coordinates of the point Q in the local system (x, y, z) can be transformed into geocentric coordinates (X , Y , Z ) using (2) (see, e.g., [3,31]): or in a shorter form using (3) where the matrices P, Q, and R are given by (4)- (6): The geocentric coordinates of P(X, Y, Z), together with an Earth Gravity Model, e.g., EGM 2008, allow for the computation of the components (ξ, η) of the deflection of the vertical (see Section 2.3 point 3 for more details). The orientation parameter Σ can be determined by solving Equation (2) for a given set of observations (s, α, β, i, j), GNSS coordinates of the station P(X, Y, Z) and the target Q(X , Y , Z ). For instance, the Levenberg-Marquardt [32] method can be used.
where the matrices P, Q, and R are given by (4)-(6): The geocentric coordinates of P(X, Y, Z), together with an Earth Gravity Model, e.g., EGM 2008, allow for the computation of the components (ξ, η) of the deflection of the vertical (see Section 2.3 point 3 for more details). The orientation parameter Σ can be determined by solving Equation (2) for a given set of observations (s, α, β, i, j), GNSS coordinates of the station P(X, Y, Z) and the target Q(X′, Y′, Z′). For instance, the Levenberg-Marquardt [32] method can be used.
The coordinates of the target A can be calculated using the set of observations obtained at 1, independently of observations obtained at 2: Since (9) and (10) should be equal, a conditional equation is expressed by (11): Equations (7), (8) and (11) constitute a system of nine independent, non-linear equations with four unknowns (s1A, s2A, Σ1, Σ2). The Levenberg-Marquardt algorithm (LMA) can be utilised to estimate unknowns. Finally, the coordinates of target A are computed from Equations (9) and (10).

Differential Relationships
A linear set of equations is required to perform the LSQ adjustment. The linear version of the equations is achieved by expanding (2) into the Taylor series. It can be achieved using the differential correction dX′ to the position X′ of the point Q (Figure 1) given by:

Differential Relationships
A linear set of equations is required to perform the LSQ adjustment. The linear version of the equations is achieved by expanding (2) into the Taylor series. It can be achieved using the differential correction dX to the position X of the point Q (Figure 1) given by: where the partial derivatives of Equations (1), (5) and (6) are given by (13)- (15):  where Finally, Equation (12) can be written as (19): The differential Equation (19) is used for setting up the linear mixed-observation Gauss-Helmert model of the spatial intersection.

LSQ Adjustment of the Intersection Based on the Gauss-Helmert Model
The non-linear LSQ solution of the intersection equations presented in Section 3.1 can be improved by considering the standard deviations of the measurements. By adding small corrections to the measured quantities as well as to the approximate values of the unknown parameters, the fundamental non-linear Equation (2) or (3) can be written in the form: Then, using Equation (19), we obtain the linear form of Equation (20): with ∆X PQ = (R P Q P P P ) T x PQ Therefore, the non-linear Equations (7), (8) and (11) can be written in linear form as (23)-(25): where v s , v α and v β are random errors of the distance and angles s, α, β, respectively, v ξ , v η are random errors of the deflection of the vertical components, ξ, η, dΣ and ds are corrections to the approximate values of the azimuth Σ and distance s, and v x = (v x , v Y , v Z ) T is the vector of random errors of the GNSS coordinates.
The system of three matrix-vector observation equations (Equations (23)-(25)) can be written as the Gauss-Helmert model in the matrix form (26): where and the vector of random errors of the observations, I is the 3 × 3 identity matrix, and p = [dΣ 1 , dΣ 2 , ds 1A , ds 2A ] T represents the corrections to approximate values of the The Gauss-Helmert model assumes condition (30): where is the weight matrix of the observations, According to the traditional Lagrange approach, the target function consisting of (30) and the linear combination (Lagrange multipliers) of condition Equation (26) is formed. The solution using Euler-Lagrange conditions (partial derivatives equal to zero) yields the following (see, e.g., [36]): where M = BP −1 B T Hence, the parameters' adjusted observations and their covariance matrices are: The covariance matrix of residual vector v is given by (37): The corrections to the observations v are acceptable if they fulfil the local test as follows (38): where i = s, α, β, X, Y, Z, ξ, η and σ v is the standard deviation of correction v. The GHM-based observation equation model is developed here to adjust the observations and parameters found in the vector p according to the LSQ method. Figure 2 shows the survey design of the experiments. The total station and the GNSS receiver were positioned at points 1 and 2 in the first experiment and at points 3 and 4 during the second experiment. Target A or targets B, C and D were located on the roof of the building some 18 m above the level of stations 1-4. The approximate distance between the instrument positions and the targets was determined using existing large-scale maps. At each station, the horizontal and vertical angles and distances were measured using a total station. The instrument height (i) and the mirror at the target were also measured. In the next step, RTK GNSS observations were taken at all points, including the target A-D points. The GNSS coordinates of the points (A-D) and the distances to them were considered as the controls of the experiments, and as such, they were not included in the calculations of the GNSS coordinates of the target points.

1.
Both experiments used a Leica Nova MS50 total station to measure the angles and distances (s, α, β). The heights of the instrument (i) and reflector (j) were also measured. The measurements obtained during the first experiment are listed in Table 1. The symbol "A" denotes an inaccessible point, while 1 and 2 denote sites of the total stations from which measurements were carried out.  Table 2.  (Figures 1 and 2a). Note that the coordinates of point A were measured as the control for the experiment.
where ϕ and λ are the latitude and longitude obtained from GNSS, and Φ and Λ are the astronomical latitude and longitude, respectively.
An alternative way to obtain the vertical deflection components is the global vertical deflection (GVD) model available from the US National Geospatial Intelligence Agency (NGA) [40]. The calculated deflection of vertical components is listed in Table 3.  Tables 4 and 5 show the observations collected while executing Experiment 2. Table 4. The total station measurement data for the testing points B, C and D (Figure 2b).  Table 5. The geocentric coordinates for Experiment 2 points B, C and D (Figure 2b). Note that the GNSS coordinates of points B, C and D were measured as controls only. Both experiments were controlled by GNSS measurements carried out at locations assumed to be inaccessible, i.e., A-D. These measurements were not a part of adjusting the observations.  (9) and (10) to calculate the coordinates of target A. Table 6 shows the results. The differences between the adjusted unknowns s(1-A) and s(2-A) and the control measurements (Table 1) equal 0.001 m and 0.021 m, respectively, while the differences between the calculated versus measured GNSS coordinates of the target A (Table 2) are −0.004 m, 0.007 m and −0.002 m. These differences are significantly smaller than three times their standard deviations.

The Gauss-Helmert Model (GHM)
The Gauss-Helmert model is used for the rigorous LSQ observations' adjustment. Unlike the LMA, it considers a priori input data accuracy, which is essential if observations are of varying precision.
We adjusted the observations using both methods to compare the LMA and GHM solutions. For the approximate values of unknowns, the LMA solution was used (Table 6). These values were also used to compute matrices A and B and the vector w according to Equations (27)- (29). Subsequently, we received the following solution for the unknown, as shown in Table 7. The GHM produced more accurate results than the LMA because the differences between the unknowns s(1-A) and s(2-A) and the controls are now 0.005 m and 0.013 m, respectively, vs. 0.001 m and 0.021 m as per the LMA solution.
This method produces more consistent results because the GHM involves all observations, coordinates of stations 1 and 2 and parameters in the adjustment. Table 7 shows point A's geocentric coordinates computed from Equations (9) and (10). The coordinates vary by 0.002 m, −0.002 m and 0.001 m from the GNSS control for X, Y and Z, respectively.
The coordinate differences of target A derived using the LMA and GHM are shown in Table 8. Evidently, the GHM method outperforms the LMA algorithm. The distance and coordinate differences are less than three times their standard deviations. The standard deviations of the computed coordinates of target A are 0.007 m, 0.013 m and 0.008 m, respectively.
The ratio of corrections of adjusted observation v and their standard deviations is 2.8, indicating that the local test-Equation (38)-has been fulfilled.

Experiment 2
Experiment 2 was conducted to further verify the GHM adjustment's usefulness in the proposed computational method for determining the precise coordinates of inaccessible locations under different GNSS satellite deployments and likely different ionosphere and troposphere conditions. The GNSS coordinates of three points, B, C and D (Figure 2b), were determined in the test.
The results show that the computed coordinates are consistent with GNSS control coordinates. The coordinate differences are below 1 cm (except for one) and are within one standard deviation of the coordinates. Also, the root-mean-square error (RMSE) calculated for the coordinate differences is below 1 cm (Table 9).

Discussion
The proposed method for determining the coordinates of inaccessible locations uses two orientation parameters (ξ, η), calculated for the stations' points from EGM 2008. Hence, the parameters link the total station measurements and the gravity field. Due to this feature, the approach allows for precise, multi-technique direct positioning in the global reference frame. The availability of the global gravity model EGM 2008 allows for straightforward calculations of (ξ, η), following Equations (39) and (40) [41]. However, to test the impact on the adjusted coordinates of the vertical deflection, we solved the case when the vertical deflection equals ξ = 0, η = 0. In this case, the rotation axis of the instrument is perpendicular to the ellipsoid. Hence, the rotation matrix Q = I, δQ/δξ = 0, δQ/δη = 0, A ξ = 0, A η = 0, meaning that the adjustment problem is simplified. We obtained target A's estimated coordinates as shown in Table 10. There is a little difference (0.002 m) in the coordinates between the simplified case (ξ = 0, η = 0) and the actual situation when, typically, ξ = 0, η = 0. The simplified solution can introduce systematic errors in the coordinates, particularly in mountainous areas or when aiming at the target requires high vertical angles. Hence, the simplified procedure is not recommended in the discussed cases.
The coordinate differences between the non-linear LMA and the GHM solutions are minor in our tests. However, the results may significantly vary in another measurement scenario (different geometry of the intersection, more/less accurate instrument) because the LMA provides a non-weighted solution compared with the GHM-weighted solution. Also, the GHM approach provides a rigorous accuracy assessment of the solution.

Conclusions
This paper introduced a method for positioning inaccessible points using the novel DIM method directly in the geocentric reference frame. This solution is based on the LSQ adjustment of measurements and parameter estimation, which includes conditional equations with unknowns within the Gauss-Helmert model. The method integrates total station observations, GNSS positioning, and vertical-axis deflection components.
The proposed method is a new contribution to precise multi-technique surveying [42]. The field tests showed that the proposed method allows for the determination of the coordinates of inaccessible locations in the geocentric reference frame at the accuracy of approximately 1 cm in terms of the standard deviation of point positioning. Many industries using these results can benefit from such accuracies. In particular, a high level of accuracy is required when monitoring the displacement and deformation of engineering structures. It can be used for engineering surveying jobs and for establishing ground control points for LiDAR and aerial photograph acquisition.
Observing the current development trends in surveying, characterised by the rapid development of modern instruments, the demand for efficient real-time or near-real-time data fusion methods is growing. These methods are vital for laser scanning and photogrammetry for object dimensioning or structural monitoring.
The proposed method is innovative, and thanks to it, one can significantly increase the accuracy and effectiveness of the job execution, mainly when a job must be performed in real or near-real time. We note such cases concerning the still-growing scope of deformation monitoring offered to construction contractors. In addition, as quantity surveying practice indicates, conducting measurement work in usually challenging construction conditions often makes it impossible to observe some control points. Such situations demand the surveyor's flexibility, and spontaneous decisions may often change the geometry of the previously planned control network in order to adapt to sudden field situations. Our proposed approach significantly reduces such risks, offering a flexible method for capturing points efficiently while keeping the highest accuracy. Future planned works include the continued development of computational methods for precise point positioning. Notably, this is important in the deformation and displacement monitoring of hazard sites near intensive industrial exploitation areas.

Data Availability Statement:
The data and the MathCad worksheet used in this study are available upon request from the co-author Edward Osada: edward.osada@dsw.edu.pl.