Design, establishment, analysis, and quality control of a high-precision reference frame in Cortes de Pallás (Spain)

A high-precision geodetic reference frame was required in Cortes de Pallás (Spain) to undertake a long-term deformation monitoring project. Involving distances up to 2000 m, the aimed accuracy was 1 mm and 3 mm for horizontal and vertical components, respectively. Such a high precision is necessary to detect possible displacements of selected points on a critical rocky area in a short period of time, i.e., 2 or 3 years, and also to provide ground control for the integration of precise 3D models periodically obtained by remote sensing techniques. Considering the historical geotechnical problems of the area, the possible influence of the hydroelectric power plant, and the peculiar orography of the zone, a proper analysis of the stability is crucial if the reference frame is to be used for rigorous over-time deformation monitoring. This paper describes the deformation monitoring of a 10-pillar geodetic network which was measured from 2018 to 2020 by using a sub-millimetric Mekometer ME5000 (0.2 mm + 0.2 ppm) along with a demanding observing methodology which includes a network of data loggers for temperature, humidity, and air pressure as well as proper calibration of sensors and instruments to prevent potential inconsistencies between the scale of the network and the unit of length of the International System (SI meter). The results demonstrated that the chosen methodology yielded the aimed accuracy. However, using such a high precision methodology entails the problem that small displacements of only 1 or 2 mm are significantly detected as deformations by conventional deformation analysis and then it arises the problem of finding a subset of stable points for a rigorous datum realization when all the points seem to displace. This general problem is analyzed in the particular case of Cortes de Pallás, where a balanced mix of deformation analysis and technical decisions was eventually adopted to define a precise and stable reference frame for rigorous over-time deformation monitoring.


Introduction
A high-precision 3D reference frame was set up as a basic component of a long-term deformation monitoring in Cortes de Pallás (Spain). To detect possible displacements of huge boulders in a short period of time, e.g., 2 or 3 years, the accuracy requirement for this reference frame was deemed to be around 1 mm and 3 mm for horizontal and vertical components, respectively. This reference frame is also used Fig. 1 Overview of the electricity power complex in Cortes de Pallás. The four extreme pillars of the geodetic network are sketched in yellow (8001,8007,8008,8010), and the critical area to be monitored is framed in red of the water pumping system (Fig. 2) could be superposed on slow terrain deformations and one-occurring motions. Four, the reference frame has to prove stable along time. Otherwise, the instability of individual points should be properly determined and considered for subsequent surveys using those points as stations.
Before pillar construction, several alternatives were analyzed for the network to meet economy, precision, and reliability in accordance with the classical approach for optimization problems in geodetic networks (Grafarend and Sansò 1985;Klein et al. 2019). Since a precise vertical control cannot be effectively included because carrying out a levelling traverse to points No.8007 and No.8010 is difficult and time consuming, the frame was carefully designed as a pure 3D trilateration network. The results obtained in years 2018, 2019, and 2020 confirmed that the aimed high accuracy was accomplished for each individual campaign.
However, a proper diagnosis of the network stability becomes crucial if the geodetic network is intended to be used as a reference frame to provide specific check points or ground control for the TDP and TLP techniques. In fact, the quality of these image-based methods as applied to deformation monitoring strongly depends on a proper Fig. 2 Detail of the La Muela II plant, with the critical area in red and four of the ten pillars in yellow integration of their 3D models or point clouds into a unique and well-defined reference system (Francioni et al. 2018).
Since optimal strategies and methods for multiepoch deformation analysis are still under discussion (Niemeier and Velsink 2019;Nowel 2019;Wiśniewski and Zienkiewicz 2016;Velsink 2015), we considered it sensible for a first diagnosis of the network stability to use conventional deformation analysis (CDA) based on the use of descriptive models (Niemeier 1981;Chen 1983;Caspary 1987;Chen et al. 1990) to test three initial hypothesis: (1) all the points are stable under a pure 3D approach, (2) only those points opposite to the critical area on the cliff (8001,8002,8008,8009) are stable under a pure 3D approach, (3) stable points are searched by an iterative 2D+1D CDA process.
Given that no conclusive results could be obtained from the three tested hypothesis, we opted for a solution that combines both CDA and geodetic criteria to determine a valid set of precise coordinates for subsequent use of the network as a reference frame as well as an assessment of the uncertainty of those coordinates.
For the sake of conciseness, only the results for hypothesis (1) and the final solution adopted will be presented and discussed.

Functional model
The measured slope distances were corrected as described in Section 2 and then adjusted by using the following functional model where, are the local coordinates of i (x j , y j , z j ) are the local coordinates of j .
The coordinates used in this functional model are referred to the CP2017 local cartesian system which was conventionally defined to facilitate the deformation analysis as well as the processing and integration of image-based techniques. The xz-plane of this local system was placed to be approximately coincident with the area of interest (Table 1).
One advantage of this conventional definition is that the coordinates and precisions obtained from the adjustments can be converted into different coordinate systems of interest with no loss of accuracy. Therefore, the adjusted coordinates for each campaign were available in the following coordinate systems: local CP2017(x, y, z), ETRS89 geodetic (ϕ, λ), and ETRS89-TM30 (E, N), along with both ellipsoidal h and orthometric H heights, thus covering all technical needs.

Individual free-network solution for each epoch
With equations of the type of Eq. 1, the sample model for epoch i is written in the usual form where, has a rank-deficiency d = 6. Therefore, the origin and orientation of the network have to be externally defined by providing c = 6 constrains (see Section "Approximate coordinates"). An alternative to overcome this datum defect is the use of the Moore-Penrose pseudoinverse (Blaha 1982), but it would entail that (1) all the points are assumed to be stable and equally valid for the datum definition, and (2) no reduction or extension of the network is admissible over time. Therefore, we opted for the use of a matrix of the datum constraints S (Blaha 1982;Caspary 1987) being (x 1 ,ȳ 1 ,z 1 , · · · ,x n ,ȳ n ,z n ) the arrangement of approximate coordinates of the n points assumed to be stable and referred to the center of gravity of the datum In order to achieve numerical stability, the matrix S t is normalized by means of S t n = (S t S) −1/2 S t such that S t n S n = I . The adjustment solution for each epoch i is then given by which is in turn used to compute the vector of estimated parametersx i , the vector of residualsv i , and the variance of unit weightσ 2 0ī Each solution was globally tested to validate the correctness of the model. Subsequently, the outlier detection was performed by using both Baarda's data snooping and Pope's Tau method (Lehman and Lösler 2016).

Procedure for testing the coordinate differences
The classical congruency test was the method used to detect possible point displacement between epochs (Niemeier 1981;Caspary 1987). This type of test aims at detecting weather or not the geometry of networks measured at different epochs can be considered statistically identical. Thus, real displacements of points can be distinguished from small displacements caused by the uncertainty of the measurements.
Assuming that two epochs are adjusted with the same datum (S-System), the null hypothesis H 0 of the congruency test is given by which is equivalent to state that the vector of their coordinate differencesd is null. The significance differencê d is tested by the well-known statistic is the covariance matrix of the difference vector h rank of the cofactor matrix Q dd σ 2 0ij =(σ 2 0i +σ 2 0j )/(f i + f j ) variance unit of weight of the combined solution.
The null hypothesis is not rejected with a level of significance α if the ratio F d follows a Fisher distribution with h and f = f i + f j degrees of freedom and the difference vectord is considered statistically null. On the contrary, if F d > F h,f,α the null hypothesis fails and it must be concluded that there has been a significant change in coordinates between epochs. The investigation of the underlying reasons for that overall network distortion is normally done by using a single point diagnosis which can be approached in different ways. We have used the conventional method that consist in splitting the quadratic form q in Eq. 9 into as many subforms as points contributing to the displacement vectord. Each single sub-form expresses the individual contribution of that point to the total form q and it is 3D tested. Alternatively, the point deformation can be tested by splitting the total 3D displacement into its horizontal (2D) and vertical (1D) components.

Site and particular aspects
The area to be monitored is a broad cliff (650-m long, 200-m high, and 300-m wide) facing the north side of a flat-topped hill in Cortes de Pallás (Spain). This critical area is located within the Cortes-La Muela hydroelectric power complex, which in turn comprises Cortes II, La Muela I y La Muela II plants.
La Muela II, near to the critical area, consists of four reversible pump turbines which use excess of energy generation capacity to store water in the upper reservoir during nighttime off-peak hours. During peak daylight hours, the accumulated water is dropped with a flow rate of 192 m 3 /s to the lower reservoir as a conventional hydroelectric facility.
Water flows between the two reservoirs via the underground power station through an 840-m long, 5.45-m diameter, and 45 • -inclined penstock housed in a conduction tunnel structure (see Fig. 2). These lifting and dropping operations are suspected to provoke vibrations and possible temporal deformations which may be superimposed on a general and slow geological deformation.
In 2015, the power plant facilities and several roads skirting the cliff were seriously damaged due to a sudden rockfall. Once the refurbishment and consolidation works were finished in late 2017, the Department of Roads and Infrastructures of the Diputació de València and the Universitat Politècnica de València (UPV) started a 3-year project to monitor both the installed anchoring systems as well as the entire cliff.

Geodetic network
A basic component of the monitoring project is a ten-pillar geodetic network (Fig. 3) whose stability is periodically monitored by using sub-millimetric EDM techniques (García-Asenjo et al. 2019). Eight of these pillars (all except No. 8007 and No. 8010 due to lack of direct visibility) are used periodically as reference stations to monitor a number of reflectors installed permanently within the critical area and to provide ground control for the TDP and TLS campaigns.
However, the area presents serious limitations for any measurement technique. Aside from the usual problems with obstacles or vegetation, there is a water reservoir which strongly limits the selection of optimal stations for both geodetic techniques and image-based techniques. The stations have to be either very close, thus impeding optimal geometries, or they have to be placed on the opposite shoreline involving distances ranging from 400 to 2000 m.

Computing process
Prior to the adjustment, the following corrections were applied to distances: refraction correction, EDM frequency drift correction, and geometric correction. Some relevant aspects are following explained.

Approximate coordinates
A set of approximate coordinates is required for both the geometrical correction and subsequent adjustment of the measured distances. However, given that the precision of the measured distances is better than 1 mm and the least squares method needs initial coordinates to be known with a relative precision of 10 −5 so that one iteration suffices, the aimed accuracy for approximate coordinates was 5 mm.
Since only distances are measured, the origin and the orientation of the 3D geodetic network have to be defined by convention. We opted for fixing the official ETRS89 coordinates of the pillar No. 8007 and the ITRF14 (2018.6299) orientation of the baselines 8007-8001 and 8007-8010, which were determined by using 10-h GNSS sessions collected during nighttime in pillars 8001, 8007, and 8010 (see Fig. 1). This initial set of approximate coordinates, which was called CP-FRM00, is also used to check and analyze all the distances for each field campaign. This pre-processing solution is carried out to detect possible problems concerning individual distances, initial coordinates or EDM scale by means of a minimum constrain solution by fixing the coordinates of points No. 8001 and No. 8007, which is the longer line, and the vertical coordinate of No. 8010. For the three pre-processing solutions, only eleven residuals showed an absolute value above 1 mm, and only one (2.71 mm) was detected as an outlier.

Refraction correction
For each distance, both the air index n and coefficient K of refraction were determined at both ends by means of the meteorological parameters. The index of refraction n was computed using the expression for the refraction index recommended in the resolution from the General Assembly of the IUGG when the required precision is 10 −7 ((IUGG 1999;Ciddor 1996;Ciddor and Hill 1999;Ciddor 2002;Pollinger 2020) and the coefficient of refraction K was obtained from the usual expression (Dodson and Zaher 1985;Baselga et al. 2014).
Once the values of n and K are determined at both ends, the average value is used to compute the first and second velocity corrections as well as the arch-to-chord correction (Bell 1992;Rüeger 1996). As expected, the last two corrections were negligible, but the absolute values obtained for the first velocity correction ranged from 5.81 to 64.74 mm.

EDM frequency correction
The ME5000 (SN 357050) was adjusted in 1989 to realize a nominal carrier frequency of 479.35817 MHz with a scale factor of 0.99999994. However, this scale factor changes with time due to the quartz oscillator ageing (Bell 1992;Rüeger 1996). Therefore, the actual frequency of the EDM has to be determined in a calibration laboratory and the derived scale factor applied subsequently to correct the field distances.
The applied scale factors were 0.918 ppm, 0.939 ppm, and 1.064 ppm for years 2018, 2019, and 2020, respectively. The rapid variation of the EDM frequency observed since 2013 might indicate a potential scale problem (Fig. 4). Unfortunately, the scale cannot be externally validated due to the lack of calibration baselines for distances of several km and nominal distances known with 10 −7 uncertainty. At present, only Nummela Standard Baseline in Finland (Jokela and Häkli 2010) could provide that standard, but its maximum length is limited to 864 m. Hopefully, the novel European standard baseline EURO5000 at the Pieniny Kippen Belt (Poland) will be available in a near future (Pollinger et al. 2021).

Geometric reductions
The measured distances are referred to the center of the instrument and the corresponding reflector. Since their height above pillars are different and vary from one station to other, distances have to be corrected to obtain the distances between the head of pillars. This geometrical correction requires the heights to be accurately measured (better than 1 mm) and good approximate coordinates. The geometric corrections ranged from −2.06 to +0.6 mm.

Solution under the hypothesis that all points are stable
Under the hypothesis that no point displaced from 2018 to 2020, the matrix of the datum constrains S t includes  (Tables 2 and 3). For all the solutions, the average values of the four corrected distances for each baseline were adjusted with a priori error of 0.2 mm + 0.2 ppm. All of them succeeded the global test of the model, being the variance of unit weight σ 2 0 = 0.439,σ 2 0 = 1.182, andσ 2 0 = 1.094, respectively. Regarding residuals, 100% were below 0.3 mm, 0.6 mm, and 0.5 mm, respectively ( Table 2).
The subsequent single point 3D analysis showed that only points No. 8001,No. 8003,and No. 8010 could be considered stable after 2 years of monitoring. A selection of the deformations obtained under this hypothesis is graphically displayed in Figs. 5 and 6.
With regard to vertical deformations, as it can be seen in Fig. 7, this hypothesis entails that the larger vertical

Solution adopted as optimal
Considering that all the hypotheses tested resulted in that most of the network points had displaced a few millimeters between the years 2018 and 2020, the underlying assumptions for the solution eventually adopted as optimal were a balanced mix of those previous results and the following technical decisions.
Firstly, we considered that the datum selection should be only based on those points that are to be subsequently used as stations to provide coordinates of check points on the rocky wall or ground control for TDP or TLS. Therefore,points No. 8007 and No. 8010 were excluded from the congruency test.
Secondly, given that none of the 3D congruency tests succeeded in finding a stable set of points, the datum selection would be specified separately for horizontal and vertical components. Consequently, the deformation analysis was also split into 2D and 1D congruency tests.
Thirdly, the criterion followed to determine the optimal datum for deformation analysis was different for horizontal and vertical components. For horizontal components, we   Regarding horizontal deformations, the CDA results confirmed that all points had significant displacements from 2018 to 2020, but only those points in the vicinity of the water reservoir reached 2D displacements around 2-3 mm (Tables 4 and 5). Moreover,points No. 8004 and No. 8005, which are in the southern shoreline, seem to be approaching points No. 8008 and No. 8009, which are placed in the northern one (Fig. 8).
Concerning vertical deformations, the decision of excluding those points seemingly sinking around the water reservoir (No. 8004,No. 8005,No. 8008,and No. 8009) proved to be valid for the definition of a stable vertical datum where significant downwards displacements can be clearly identified (Fig. 9). The subsequent single point 1D analysis showed that all points defining the vertical datum can be considered stable from 2018 to 2020 (Tables 6  and 7).

Conclusions
In the light of our experience in Cortes de Pallás, some conclusions can be drawn. First, geodetic techniques, and specifically sub-millimetric EDM based ones, can provide coordinates with a precision better than 1 ppm even at distances larger than 1 km. Therefore, the establishment of well-controlled reference frames in critical deformation monitoring projects is not only possible, but also mandatory. Otherwise, long-term monitoring projects based on periodical campaigns by using traditional surveying or remote sensing techniques cannot be consistently integrated in a unique datum and thus, the deformation diagnosis derived from them will be consequently biased.
Second, when sub-millimetric geodetic techniques are used, especially in areas with apparent geotechnical problems like the current case in Cortes de Pallás, the comforting hypothesis of having stable points to serve as reference frame is no longer valid because small horizontal displacements of only several millimeters can be significantly detected in a short period of time with the mere use of the conventional deformation analysis. Consequently, the question of finding a stable datum becomes mathematically problematic and no clear procedure can be found in literature. In this case, we opted for a solution suitable to the reference frame practical purpose, although it is clear that it is a matter of time that the approximate coordinates obtained in 2018 do not fulfil the requisites for conventional deformation analysis and new proposals should be investigated. Third, although it is commonly admitted that traditional surveys based on total stations or 3D point clouds derived from remote sensing techniques are reasonably well integrated by using GNSS techniques, our experience demonstrates that monitoring of high-precision reference frames requires state-of-the-art geodetic instruments, like long-range sub-millimetric EDMs, as well as length metrology methods.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.