Development of a Submillimetric GNSS-Based Distance Meter for Length Metrology †

Absolute distance determination in the open air with an uncertainty of a few tenths of a millimetre is increasingly required in many applications that involve high precision geodetic metrology. No matter the technique used to measure, the resulting distances must be proven consistent with the unit of length (SI-metre) as realized in the outdoor facilities traditionally used in length metrology, which are also known as calibration baselines of reference. The current calibration baselines of reference have distances in the range of 10 to 1000 m, but at present there is no solution on the market to provide distances with submillimetric precision in that range. Consequently, new techniques such as multi-wave interferometry, two-wave laser telemeters or laser trackers are being developed. A possible alternative to those sophisticated and expensive techniques is the use of widely used Global Navigation Satellite Systems (GNSS) in order to provide a GNSS-Based Distance Meter (GBDM). The use of a GBDM as a potential technique for length metrology has been thoroughly analysed in several European research projects by using the state-of-the-art geodetic software, such as Bernese 5.2, but no definite conclusions have been drawn and some metrological questions are considered still open. In this paper, we describe a dedicated approach to build up a submillimetric GBDM able to be applied in the current calibration baselines of reference, as well as possible methods to cope with the multipath error of the GNSS signals which is the major limitation for the practical uptaking of the technique in metrology. The accuracy of the proposed approach has been tested following the length metrology standards in four experiments carried out in the Universitat Politècnica de València (UPV). The results demonstrate that the proposed GBDM can provide an accuracy of a few tenths of a millimetre in the current calibration baselines of reference.


Introduction
Absolute distance determination in the open air with an uncertainty of a few tenths of a millimetre is increasingly required in fields such as length metrology, high precision geodetic metrology, deformation monitoring at critical sites or local ties in Fundamental Geodetic Observatories (FGOs) [1][2][3].
In the present context, the term absolute refers to the fact that distances have to be consistent with the metre, the unit of length of the International System of Units (SI) which is currently defined as the length of the path travelled by light in vacuum during a time interval of 1/299792458 of a second [4]. This definition of SI-metre is realized by the National Metrology Institutes (NMIs) with a standard uncertainty U (k=1) = 10 −12 by using laser interferometers with stabilised frequencies, though such a high accuracy can only be obtained for short distances, i.e., several metres, under laboratory conditions. Following the core concept of metrological traceability, those references or primary standards are subsequently transferred or disseminated from laboratories to outdoor metrological facilities that are also known as calibrations baselines. These calibration baselines, whose inter-pillar lengths usually range from 10 m to 1000 m, are then used as secondary and working standards for practical purposes [5][6][7].
The only calibration baseline that has proved stable with a relative standard uncertainty near u (k=1) = 10 −8 for more than 70 years is the Nummela Standard Baseline (NSB) in Finland (864 m), and thus acknowledged as international standard in geodetic metrology. Such a high accuracy can only be obtained by using the Väisäla interferometer which has the additional metrologic advantage of being directly traced to SI-metre [8][9][10]. With this method being very demanding, the NSB inter-pillar distances are measured only every five years by the Finnish Geospatial Research Institute (FGI), although its accurate scale can be traceably transferred to other baselines of reference by using submillimetric Electronic Distance Meters (EDM) like the Kern Mekometer ME5000 (ME5000) which was manufactured in the late 1980s [11,12].
Nevertheless, the use of the out-of-production ME5000 entails two problems. Firstly, only a few units of this type of EDM are still in use, and secondly, having only one carrier phase (0.628 nm), the air refractivity has to be determined with an accuracy better than 0.1 ppm by using meteorological sensors, which is rather difficult for distances longer than 1 km [12]. Consequently, new techniques such as multi-wave interferometry, two-wave laser telemeters or laser trackers are currently required to go beyond the ME5000 [13][14][15]. Yet the research effort is considerable, and it demands joint projects like the current 18SIB01 GeoMetre-large-scale dimensional measurements for geodesy' which has received funding from the EMPIR programme co-financed by the participating states and from the European Union's Horizon 2020 research and innovation programme [16].
A possible alternative to those expensive techniques is the use of Global Navigation Satellite Systems (GNSS) to develop a GNSS-Based Distance Meter (GBDM). Some experiments that compare GNSS distances with EDM distances can be found in the literature, but being most of them based in the use of commercial GNSS software and EDM instruments, the differences between these two techniques can easily reach around 3 mm [17,18]. The GBDM as a potential technique for length metrology was firstly analysed in the EMRP-JRP 3.1-TP 3 'Length' [19], and subsequently in the EMRP-JRP SIB60 'Surveying' [20]. Most of the studies were performed by using state-of-the-art geodetic software, i.e., the Bernese 5.2 [21], along with the most advanced geodetic products and models. The reporting document for the GBDM technique recommends the use of carrier-phase double differences (DD), and points out that GNSS signal multipath is the major limiting factor to get submillimetric accuracy, but no recommendation is given concerning a possible standard strategy for optimal GNSS processing and some metrological questions such as uncertainty determination or traceability are considered still open [22]. Consequently, some criticism still remains for the GNSS to be used as a submillimetric technique for length metrology.
However, previous research conducted at the Department of Cartographic Engineering, Geodesy and Photogrammetry of the Universitat Politècnica de València (DICGF-UPV) showed the potential of GNSS techniques in terms of reproducibility for measuring distances over several hundred metres with uncertainties below one millimetre by using a functional model different from the traditional geodetic approach along with robust estimation [1, 23,24]. The initial assumptions, the functional models and the intended computing approach are described in Section 2.
As previously mentioned, the GNSS signal multipath error is considered the major limitation to use a GBDM in current metrologic facilities. The traditional geodetic approach to cope with this problem is to carefully select the GNSS sites and average out the multipath error by using long observation span times. However, current calibration baselines of reference [10,25] were mostly set up decades ago, and they do not always have favourable conditions for GNSS measurements since they are normally surrounded by trees or buildings. As a result, L1/E1 DD residuals can easily reach up to several centimetres with no Gaussian distribution and consequently the obtained distance will be biased, thus impeding the desired submillimetric accuracy. Therefore, a consistent GBDM should include some processing strategy for multipath mitigation. Therefore, the last field experiment was specially designed to investigate and develop possible methods for multipath mitigation as applied to length metrology. This part of the investigation was conducted at the Space Geodesy and Navigation Laboratory of the University College London (SGNL-UCL), and their results are shown in Section 3.
Nonetheless, not only does the GBDM technique have to be proven to be precise but it must also be accurate according to length metrology standards. Thus, the proposed GBDM technique was tested in four field experiments conducted at the UPV submillimetric calibration baseline from the years 2013 to 2016. Each experiment was a four-day campaign where both the GNSS and EDM techniques were applied. During the daylight hours, the inter-pillar distances of the calibration baseline were measured using a ME5000 in accordance to ISO17123-4 [5], and subsequently, GNSS antennas were set up to collect twelve-hour sessions overnight. The same operational scheme was repeated for four consecutive days. The experiment was designed in such a way that the obtained distances can be traced to the SI-metre.
This paper is organised as follows. In Section 2, the general approach for GNSS data processing is described. Section 3 is devoted to the methods for multipath detection and mitigation that were developed and tested. The assessment of the resulting accuracy by comparison with distances derived from EDM traced to the SI-metre definition in a calibration baseline is explained in Section 4. Finally, some conclusions and future work are summarised in Section 5.

Basic Carrier-Phase Equation
Unlike the traditional EDM distance determination, where electromagnetic signal path travels several hundred metres between both ends of the measured baseline, the GNSSbased distance has to be indirectly obtained from electronic signals that travel from satellites to the antennas which involves distances around 22,000 km. The general equation for a carrier phase observable can be written as follows where: Φ k i (t) carrier phase measurement collected at receiver time t i , in m ρ k i (t) geometric distance between the receiver i and the satellite k, in m c speed of light in vacuum, in m/ṡ ρ k i (t) topocentric range rate of the satellite k with regard to the receiver i, in m/s dt i receiver clock offset, in s dt k satellite clock offset, in s carrier phase wavelength, in m N k i integer ambiguity for the carrier phase observable, in cycles Φ i (t 0 ) initial carrier phase offset for the receiver i, in m Φ k (t 0 ) initial carrier phase offset for the satellite k, in m δ i receiver hardware delay, in m δ k satellite hardware delay, in m δm k i (t) carrier phase multipath error, in m k i (t) random noise, in cycles For the sake of simplicity, the notation assumes the use of only one constellation (e.g., GPS).

Initial Assumptions
Our approach to investigate the capabilities of GNSS for length determination with submillimetric accuracy relies on four initial assumptions: First, modern receivers provide carrier phase observations with submillimetre level noise. Secondly, distances determined by GNSS are highly stable at a global scale. Thirdly, the modulus of a vector between two GNSS receiver antennas can be more precisely determined than its corresponding components. Finally, the increasing number of operational GNSS satellites allows one to optimise the selection of satellites to improve length determination. These four assumptions are discussed below.
Regarding the noise of the GNSS carrier phase signal, it used to be assumed as a random error of 1% of the carrier wavelength, but recent literature shows that modern GNSS receivers can provide enhanced performance [26,27]. The random noise of a GNSS signal depends on the electronics of the receivers, but also on the quality of the signal. Therefore, a different noise should be expected for each constellation and even for each individual satellite. A standard method to assess the overall random noise of a GNSS receiver is to perform a zero baseline test, which consists of processing the double differences of phase observations of two receivers that are connected to the same antenna by means of a signal splitter (see Figure 1).  Figure 2 shows the double difference residuals obtained using two multi-constellation Leica GS10 receivers. The overall root mean square (RMS) value obtained for L1 double difference residuals was 0.7 mm, thus confirming our first initial assumption that they are currently below 1 mm × 2, that is, as we assumed, less than 1 mm of standard error for a carrier phase observation amplified by a factor of 2 due to the double difference combination (see e.g., expression (20.67) in Reference [28]). Concerning the GNSS scale stability, it is maintained worldwide at the level of 1 ppb = 10 −9 , which is to say 1 mm over 1000 km. Nonetheless, some open questions still remain as to whether this high absolute scale accuracy is also applicable to short distance determination in the existing metrological infrastructures. Moreover, the GNSS scale stability relies on the accuracy of the precise ephemeris and ultimately on the International Terrestrial Reference Frame (ITRF) which is used for their determination. Since the ITRF scale is determined by a combination of Satellite Laser Ranging (SLR) and Very Long Baseline Interferometry (VLBI) observations whose scale differs by about 10 −9 [29], the linkage between the centre of the instruments used by SLR, VLBI and GNSS at Fundamental Geodetic Observatories (FGO), also known as local ties, should be performed with submillimetric accuracy [2]. Thus, an additional benefit of submillimetric length determination with an appropriate metrological uncertainty evaluation based on traceable instruments and methods would be to improve the ITRF scale and consequently the accuracy of the precise ephemeris of GNSS satellites.
The third assumption is that GNSS distance may be better determined than the individual components of the corresponding vector. Double difference processing usually yields precise vectors whose components ∆X, ∆Y, ∆Z are subsequently used to compute the distance Fortunately, correlation usually leads to σ D < σ 2 ∆X + σ 2 ∆Y + σ 2 ∆Z (and even σ D < σ ∆X and σ D < σ ∆Y , σ D < σ ∆Z ) because in the expression the covariances may be positive or negative depending on the relative geometry between the baseline and the observed satellites. Thus, a proper selection of observables may lead to a reduction in the standard error of the obtained GNSS distance. Finally, a multi-constellation approach has clear advantages for length determination if compared to a single-constellation one, i.e., only GPS. A large number of satellites in view allows the inclusion in the system of normal equations of only those double differences whose impact on the baseline is lower than an assumed threshold as well as to eliminate those satellites whose phase observables contain a strong multipath error. As an illustration, in the last GNSS campaign that was carried out at the UPV calibration baseline in July 2015, the number of satellites registered in a time span of eight hours was the following: 21 GPS, 19 GLONASS and 4 Galileo. Thus, the inclusion of the new and emerging constellations gives rise to new possibilities for precise GNSS length determination.

Functional Model for the Initial Coordinates
The GBDM technique starts with the determination of the absolute coordinates of the ends of the measured distance. This initial solution has to provide three-dimensional coordinates with an accuracy of a few centimetres in the same coordinate reference frame as the precise ephemerides used in the second step. This requirement is not very difficult to fulfil, and can be achieved by using different strategies, but we recommend the use of a common Precise Point Position (PPP) solution. This approach, which is based on an epoch-wise processing strategy that uses an ionosphere-free combination with floating ambiguities, is relatively easy to implement or can be alternatively obtained from online services like CSRS-PPP [30].
Some advantages of this straightforward approach are that the bulk of the ionosphere error is removed, slow changes of the tropospheric delays can be easily absorbed by parametrization, the receiver clock offsets are well determined and the obtained ionospherefree residuals are expected to be mostly multipath error and possible antenna calibration mismodelling.
Concerning the functional model, a widely used equation for receiver i and satellite k is (see, e.g., [31]) where: elevation and azimuth angles (rad) T zwd zenith wet troposphere delay (m) m G () Global Mapping Function (GMF) G N , G E north and east troposphere gradients The non-integer ambiguity term A k i comprises the integer ambiguity, the initial phase offsets for both frequencies, and the hardware delays. The vector of Kalman filter states that are estimated is Taking into account the implemented PPP approach, the obtained ionosphere-free residuals are expected to be mostly multipath error and possible antenna calibration mismodelling [32].

Functional Model for the Double Differences
The functional model used in our approach is specifically tailored to length determination and differs from the corresponding geodetic approach in four ways: only L1/E1 double differences are processed, ambiguity determination is not required, and both ionospheric and tropospheric corrections are not included.
Since L1/E1 carrier phase observables are less noisy than L2, L5, E5a or E5b observables, or any combination of them [33], double differences of L1/E1 carrier phase are the preferred equations to obtain the final solution.
A free ambiguity approach reduces the number of unknowns to determine thus strengthening the model [34]. The only requirement for using the free-ambiguity method is a prior determination of the absolute coordinates of both ends with an accuracy of a few centimetres, which can be achieved by means of a PPP static processing, thus obtaining (X for substraction with the order and expanding around (X i ) held fixed, Equation (6) can be rewritten as which can be subsequently arranged as where the only unknowns to be solved are dX j , dY j , dZ j . Regarding tropospheric delays, our current model assumes that they are negligible for distances below 1 km. Note that existing metrological facilities are small and points of interest are usually at about the same altitude. Consequently, the path taken by a GNSS signal through the atmosphere to each end of a baseline is almost the same, especially for those satellites with higher elevation. In addition, we collected GNSS data in the nighttime. Thus, only residual atmospheric delay errors with minor influence on the distance are expected.
For longer distances, ionospheric delay can be determined using a combination of carrier waves [28] and subsequently use them to correct L1/E1 carrier phase observables. On the contrary, the influence of residual tropospheric delays should be taken into account when both ends of the measured baseline have a substantial height difference which can happen in some local ties or geodetic control networks. Since the tropospheric errors expected in the double difference equations are very small and also dependant on local meteorological parameters, global models may not be accurate enough for submillimetric length determination.

GNSS Multipath Mitigation
The GNSS multipath error arises from the fact that the signal from a satellite arrives at the antenna via multiple paths due to reflection and diffraction. These indirect-path signals can distort the received signal by up to one quarter-cycle and this error does not cancel or mitigate when forming the double difference combination. Therefore, GNSS multipath becomes a critical source of error in high precision applications [28,35], and even more for length metrology because most of the current calibration baselines of reference, which were set up decades ago before the advent of the space geodesy era, are surrounded by objects such as trees, buildings or other facilities.
Multipath interference can be separated into near-field and far-field effects, which do have different properties. Near-field multipath results from the closest vicinity of the antenna, i.e., the first 50 cm around the antenna, and mainly leads to a systematic bias that can be prevented by creating identical near-field situation in all stations which involves the same type of pillars, GNSS antenna model, antenna mounting, and routing of the antenna cable. However, the near-field can also change the electronic properties of the antenna and thus, individual antenna calibrations can be considered accurate as long as the actual near-field situation was reproduced during the calibration procedure. Nonetheless, this refinement is not usually included in the calibration process due to size and weight limitations and it would only make sense if each calibrated antenna was always setup in the same pillar [22,[35][36][37].
On the other hand, far-field effects arise from reflecting surfaces in the environment of the antenna such as trees or buildings, and for static surveys, which is the case at hand, lead to short periodic errors. Following the traditional geodetic approach, the far-field multipath is generally assumed to average out by sufficiently long observation times [35,38,39]. In consequence, this approach is particularly valid for permanent GNSS stations which have a favourable observational environment and 24-h sessions but cannot always be applied to calibration baselines and short time spans. Our initial approach to cope with multipath was based on the use of the signal-to-noise ratio flag and global robust estimation [40].
Interestingly for our study, the trees surrounding the UPV calibration baseline, which were young in 2012, are increasing in size as the years go by. As a consequence, our former approach to cope with multipath, based on the use of the signal-to-noise ratio flag and global robust estimation [40], cannot be successfully applied under these new conditions. To give a picture of the problem, some DD residuals which were obtained in our experiment are plotted in Figure 3. It can be clearly seen that the L1 DD residuals range from 0.04 m to −0.05 m and they cannot be considered Gaussian. Therefore, the obtained distance will be biased. In addition, the L2 DD residuals are bigger, as expected, and follow a different pattern than the L1 DD residuals. Similar residuals were found in other DD combinations.
Over the last 10 years, the Space Geodesy and Navigation Laboratory of University College London (UCL) have developed and successfully applied techniques for GNSS multipath detection and mitigation, some of them suitable for GNSS processing for distance determination [31,32,41,42]. Specifically, we have now implemented the observationdomain sidereal filtering technique and the signal-to-noise ratio (SNR) estimator.
The sidereal filtering technique [31,32] takes advantage of the consecutive passes of each satellite to predict the multipath error. As long as both azimuth and elevation repeat for the same satellite, the multipath effect can be considered similar, and subsequently corrected. For GPS satellites, the repeating time is one sidereal day, which makes the technique attractive to use if data from consecutive days are available, whereas for Galileo satellites the repeating time increases to 10 sidereal days and the use of the technique would require GNSS observation sessions 10 days apart, which is impractical. In what follows, we will restrict the application of the sidereal filtering technique to GPS observations only.
In contrast, the SNR estimator does not provide a value for the multipath error [41]. By differencing SNR measurements on different frequencies and comparing the result with that obtained in a low-multipath environment, multipath can be detected. The result is a number that can be used, in theory, to identify the presence of multipath. To characterize a low-multipath environment, we opted for the VALE permanent station. It is only 100 m away from the UPV calibration baseline and its equipment (Leica GR10 receiver and Leica AR25.R3 antenna) is very similar to the equipment used in our experiment described in the next section. Figure 4 shows the SNR estimator for the Galileo satellite that was selected as a low-multipath reference. The line in black is the expected value, and in green, orange and red are respectively the 1σ, 2σ and 3σ thresholds. If the estimator is clearly above the red line, multipath is considered to be present. Both techniques, sidereal filtering and SNR estimator, may be complementary. This is important since each technique may on certain occasions yield null values even when multipath is present. Conversely, depending on some aspects of the PPP implementation, e.g., cycle-slip detection and correction strategies, the residuals obtained could be abnormally high with no relation to multipath. However, matching residuals by using azimuth and elevation has clear advantages over time. First, the technique is independent of the value of multipath and does not need to rely on the identification of the correct satellite repeat period (a significant problem so far, as pointed out in Reference [43]). Second, the multipath model could be eventually parameterised in terms of azimuth and elevation even for the satellites with low multipath.
In summary, the whole process to mitigate multipath in the context of GBDM consists of the following four steps 1. To obtain a PPP solution for each site and for every day using an ionosphere-free combination. 2. To match the ionosphere-free observation residuals on successive days using azimuth and elevation in order to obtain a site-specific and time-specific multipath model for an ionosphere-free combination. 3. To form ionosphere-free DDs using the obtained multipath models with the same DD scheme that will be eventually used in the final L1/E1 solution. 4. To determine the L1/E1 DD multipath corrections from the previous values by using DD residuals from an independent DD solution and the same ionosphere-free combination that was used to obtain the PPP solution.
In the first step, the PPP solution is obtained with an epoch-wise processing strategy that uses an ionosphere-free combination with floating ambiguities. Some advantages of this straightforward approach are that the bulk of the ionosphere error is removed, the multipath error is amplified with respect to that of the original observed frequencies, and it can be implemented in real-time if desired.
Nevertheless, part of the multipath error could slip into other parameters and errors caused by the imperfect modelling of other parameters (receiver clock, troposphere delay, phase ambiguities) could be wrongly considered as multipath. For that reason, it is relevant to look for a repeating pattern in the second step. Only the residuals that follow a sidereal pattern can be considered as multipath error. Once the consecutive days are matched using azimuth and elevation, the median of the residuals is retained as multipath error. The result is a site-specific and time-specific multipath model for the combination used in the PPP solution. Interestingly, the model formed may contain some possible antenna calibration mismodelling since this error also repeats with azimuth and elevation.
In the third step, the site-specific and time-specific models that were built up for each satellite are combined to form the required ionosphere-free DD residuals. This combination largely cancels those components of the residuals which originate in the satellites, such as orbits and clocks, those related to the receiver clock offsets and possible residual atmospheric errors.
Finally, the L1/E1 DD multipath corrections are obtained for each day by using the ionosphere-free DD corrections and the corresponding DD residuals from an independent DD auxiliary solution. For that independent DD solution several frequencies such as L2, L5, E5a or E5b can be used as long as they were employed in the PPP solution. Obviously, the noise of the new multipath-corrected DD equation is higher then the original one, but that is not a problem for the case at hand as long as the multipath error is mitigated in such a way that the resulting distance is not biased.
The observation-domain sidereal filtering for GPS Precise Point Positioning (PPP) developed at the SGNL-UCL [31] uses the ionosphere-free linear combination L 3 with the following combination factors Please note we have chosen the notation L 3 , which is used in the Bernese software, to be consistent with L 1 and L 2 .
The formed L 3 combination has a wavelength λ 3 = λ 1 λ 2 77λ 2 − 60λ 1 0.006 m and with standard noise values σ L 1 = σ L 2 = 2 mm, the resulting L 3 noise is Nonetheless, modern receivers can measure the carrier phase with better precision and from our experience the values σ L 1 = σ L 2 = 1 mm are more realistic so that the resulting noise for the L 3 combination can be considered around 3 mm.
Since the noise of the ionosphere-free combination is σ L3 ≈ 3 mm, the obtained DD residuals are expected to have a theoretical noise σ DD L3 ≈ 6 mm that could hamper the determination of a L 1 site-specific multipath model by using a PPP approach. Surprisingly, if DD L1 residuals are obtained straight from the DD L3 residuals by substracting the corresponding DD L2 residuals, the resulting residuals have a noise σ DD L1 ≈ 2 mm. The reason is that the coefficients are smaller than one.
and therefore, Figure 5 illustrates the PPP ionosphere-free phase residuals for pillar 3 (PL3A) (see Figure 1) and satellite G12 after the matching process. As can be seen, the residuals from the four consecutive days of the experiment are consistent and show a clear sidereal pattern. Therefore, the initial assumption for the PPP residuals to be mostly multipath error is confirmed. When the residuals from the consecutive days are compared with each other, the differences are normally below 5 mm, thus confirming the correctness of the functional model used in the PPP implementation and its suitability for this application. Nevertheless, in some cases the PPP residuals differ by up to 1 cm, possibly due to a slightly different satellite trajectory, different humidity of the ground or technicalities of the particular PPP implementation. In those cases, there is still some consistency between days, and normally only one day at a maximum is different from the rest, so that retaining the median as the multipath value for the model usually solves the problem. The figure above is a zoomed-in excerpt from Figure 6, which shows the residuals and the multipath model for the entire observation time for this particular satellite, that is, the approximately 6.5 h the satellite has been with an elevation above 10º. Evidently, the observations of low elevation are strongly affected by multipath and should not be used for precise positioning. Now let us focus on the impact on DD processing. In Figures 7-9 the difference between the initial L1 DD residuals obtained from the Bernese solution (blue) and the L1 DD multipath correction given by the proposed multipath mitigation method (green) are plotted (red) for three different DD combinations and days. In the three examples, the L1 DD multipath error is mitigated and the new observed minus computed value is lower than one centimetre on average. The three examples also show that the proposed method seems to work independently of the elevation of the satellites, and even in the presence of data gaps.   As it can be seen Figure 7, the range of the obtained differences is 5 mm on average. However, some values at around 20.6 h have an absolute value of about 2 cm with residuals peaking −5 cm. Even in that case, the influence of the multipath error is largely diminished in the new observed-minus-computed-term, and the overall impact on the distance obtained can still be considered unbiased.
The example in Figure 8 shows higher residuals that reach up to 8 cm. In this case, the bandwidth of the differences is over 1 cm, and slightly higher if compared to the previous example. Again, the differences seem to be randomly distributed around zero and no impact is expected in the subsequent obtained distance.
Finally, the example given in Figure 9 has been chosen because the differences between the initial L1 DD residuals and the multipath corrections seem to be not so well centred as other DD combinations. In this case, although the overall range of the differences is over 8 mm, it can be seen that it is wider and slightly biased when both satellites are at a low elevation.

The UPV Calibration Baseline of Reference
To validate the accuracy of the GNSS-based distances they have to be compared with distances accurately traced to SI-metre in an outdoor metrological infrastructure. For practical reasons, we opted for the UPV calibration baseline. This baseline, which was set up in November 2007 by the Department of Cartographic Engineering, Geodesy and Photogrammetry (DICGF), was originally planned to evaluate the uncertainty of geodetic instruments according to ISO 17123 series [44]. In 2012, absolute scale was transferred from de Nummela Standard Baseline (NSB) by the Finnish Geospatial Research Institute (FGI) so as to carry out the experiment to validate the accuracy of the GNSS-based distances.
As shown in Figure 10, the UPV metrological infrastructure is a triangle-shaped test field with seven observation pillars which includes an Heerbrugg-type EDM calibration baseline [12].  The pillars have a diameter of 22 cm and height 1.20 m above ground level made of two insulated steel pipes. The inner one covers a concrete structure with an square metre foundation that extends to a depth of 60 cm. The outer steel pipe prevents the inner pillar from differential dilations due to meteorological effects. To install measuring instruments on the top they have a double forced-centring mount system: the standard 5/8" fixing screw and a Kern-type trivet system. Three measurement campaigns were carried out from December 2007 to February 2008 using the following total stations: Topcon GTS605, Leica TPS 1200+ and Leica TDA5005. All of them proved to be compatible within a combined least squares adjustment and the coordinates given in Table 1 were adopted as first set of approximate local coordinates. The coordinates are referred to a local right handed Cartesian system whose origin is the top of the pillar No.1, the x-axis points to pillar No.6 and z-axis is positive upwards. The horizontal and vertical stability of the pillars over time has been periodically monitored [45][46][47]. Since no significative displacement was detected, this first set of coordinates have been safely regarded for geometrical reduction computations in subsequent campaigns.
The field campaign to transfer scale from the Nummela Standard Baseline (NSB) in 2012 was carried out from 28 May to 1 June by the Finnish Geospatial Research Institute (FGI). The measurements were done using the Kern ME5000 No. 357094 of the Aalto University along with the reflector No.374414 which were calibrated at Nummela before and after the transferring campaign. Dry and wet temperatures were measured using calibrated Thies Clima Assmann-Type psycrhometers with an estimated ±0.3 • C uncertainty and air pressure was measured using two calibrated Thommen 3B4.01.1 aneroid instruments with an estimated ±0.3 hPa uncertainty. The EDM and the reflector were installed in the 5/8" fixing screws of the pillars using Leica GDF321 tribrachs. As the observational scheme consisted in four sets of 'double-in all-combinations'distances, every distance was consequently measured 16 times. This observational scheme was strictly repeated in the subsequent four field campaigns.
Therefore, the traceability chain for the UPV baseline distances that were to be use as ground truth was: SI definition of the metre, quartz gauge bar, Väisälä interferometer, FGI Nummela Standard Baseline, FGI Mekometer ME5000 and UPV calibration baseline. The resulting nominal distances were obtained with expanded total uncertainties U k=2 of 0.1-0.3 mm (Table 2).

Description of the Experiment
To be prove accurate, the GNSS-based distances had to be compared with the SItraced distances provided by the FGI in 2012. Furthermore, since accuracy comprises both repeatability and reproducibility, the comparison had to be done over time with GNSS measurements carried out under different environmental conditions. Nonetheless, the SI-traced distances provided by the FGI were realized by the pillars in a specific epoch (June 2012), and consequently their stability over time has to be considered as a critical part of the study. Considering the construction features of the pillars, errors due to vibrations and reproducibility of the centre mount can be assumed negligible, but the possible instability of pillars, which chiefly depends on the grounding and the geological profile, would diminish the uncertainty of SI-traced distances of reference. The instability of pillars is normally included in the corresponding uncertainty budget [22], but detecting possible pillar displacements at submillimetric level requires special instruments like the out-ofproduction ME5000 or laser trackers [14,48]. For this reason, the experiment was planned to be coincident with the deformation monitoring campaigns of the UPV calibration baseline which were annually carried out from year 2013 to year 2016. In each campaign, the EDM and GNSS techniques were interlaced in a way that the distances obtained can be assumed to be done at the same epoch. The EDM distances were measured in daytime, while the GNSS observations were collected in 12-hour sessions during nighttime. In the following paragraphs, the operational scheme of both types of measurements are further described.
With regard to the EDM campaigns, an agreement between the UPV and the Universidad Complutense de Madrid (UCM) enabled us to use the Mekometer ME5000 No. 357094 along with reflectors No.374447 and No.374448 to carry out the deformation monitoring of the UPV baseline [46]. The dates and instruments used in each campaign are shown in Table 3. Dry and wet temperatures were always measured at both ends by using calibrated Thies 4410 Assmann-Type psychrometers with an estimated ±0.3 • C uncertainty and air pressure was measured using two calibrated Thommen 2A4.611 aneroid instruments with an estimated ±0.3 hPa uncertainty. The EDM and the reflector were installed in the 5/8" fixing screws of the pillars using Leica GDF321 tribrachs. Each campaign consisted in four series following exactly the same observational scheme as the first campaign carried out by the FGI in 2012 to transfer absolute scale from Nummela. Since each series of EDM measurements takes around one day, the total observation time for each campaign is four days. This observational scheme not only provides the actual EDM distances between pillars with uncertainties of 0.1-0.3 mm, but also the zero-error for the set distance meter and reflector, which is determined in accordance to the full procedure of the ISO17123-4. The raw distances measured in the four field campaigns were processed following exactly the same steps that the FGI scale-transferring campaign in 2012. The three corrections that were sequentially applied were: refraction correction, calibration correction and geometrical reduction.
The information required for computing the refraction correction consists of dry temperature, wet temperature and air pressure collected at both ends of the beam path for every measured distance along with the calibration values of the thermometers and barometers used in the field campaign. The actual refraction index for every single distance is computed using the precise Ciddor and Hill method [49,50] as recommended by the International Association of Geodesy when the scale is required to be better than 1 ppm. It is worth noting, however, that this formulation has an incorrect sign, as recently pointed out by [51]. Once obtained the actual refraction index, the three meteorological corrections recommended for the ME5000 [11] were computed, though only the first speed correction was retained because the second speed correction and the beam correction were negligible even for the longest distance (330 m).
The calibration correction was then applied to the slope distances resulting from the refraction correction. Applying the calibration parameters (zero-error and scale error) improves the accuracy of the measured distances, though the precision of the resulting distances decreases as the uncertainty of the calibration parameters has to be taken into account.
Finally, the geometrical reduction is applied to obtain the distances in Table 4. Concerning the GNSS campaigns, they were carried out in 12-h static sessions during the nighttime once the daily EDM measurements were finished. Table 5 shows the dates and other details concerning the GNSS campaigns. Since pillars No.1 and No.3 had better GNSS observing conditions than the rest, they were selected to conduct the experiment. Nonetheless, the trees surrounding the baseline had considerably grown from 2013 to 2016, making the multipath error apparent even for the selected pillars. Therefore, we decided to reduce the GNSS data interval from 15 s to 1 s in the 2016 campaign so that the methods to detect and mitigate the multipath error that were being developed at Space Geodesy and Navigation Laboratory of the University College London (SGNL-UCL), which were described in the previous section, could be better implemented and tested. With regard to the GNSS antennas, we always used the same two 3D choke-ring Leica AR25.R4, which were individually calibrated by the University of Bonn by using the anechoic chamber method [52,53]. The two antennas were set up carefully oriented towards the geodetic north which was rigorously realized by benchmarks that are permanently located near every pillar. The GNSS distances finally obtained are given in Table 6. It is worth noting that all the DD residuals and testing results in this study are consistent with the solution by Bernese 5.2 [21]. The corresponding differences among the GBDM distances, the ME5000 distances, as well as the original calibrated distance traced to the metre-SI (which is displayed in Table 2) are shown in Table 7. The GNSS measurements in 2013 were collected using four single-constellation (GPSonly) Trimble 5700 receivers and two signal splitters. They were performed in the daytime shortly before the corresponding monitoring campaign using the UCM Mekometer ME5000. The differences between the computed GNSS distances and their corresponding FGIcertified and UCM-Mekometer ME5000 were respectively +0.12 mm and −0.13 mm. Taking into account their respective uncertainty (standard deviation for GNSS), the difference can be considered statistically negligible. In addition, there is no significant evidence of pillar displacements from 2012 to 2013.
In 2014, the GNSS measurements were collected using two single-constellation Leica GS10 receivers and they were performed in the nighttime whereas the corresponding monitoring campaign using the UCM Mekometer ME5000 was carried out during the same days in the daylight. The comparison between the GNSS distance and the EDMderived for this epoch yielded a difference of 0.28 mm, which can be considered statistically negligible taking into account their respective uncertainties (standard deviation for GNSS). In addition, non-significant is the baseline length variation in this epoch with respect to the previous epoch. This agrees with the conclusion obtained by a dedicated deformation approach for the entire network [47], which concluded no significant displacements for the baseline.
The GNSS measurements in 2015 were collected in the nighttime using two multiconstellation Leica GS10 receivers. The spantime was eight hours and the receivers measured carrier phase observables for the following number of satellites: 21 GPS, 19 GLONASS and 4 Galileo. The corresponding monitoring campaign using the UCM Mekometer ME5000 was carried out during the same days in the daylight. Since both EDM-based distances, i.e., FGI-certified and the one obtained by using the UCM Mekometer ME5000, agreed well, it is also concluded that there were no significant displacements. Once again, the difference between the obtained GNSS distance and the corresponding using the UCM Mekometer ME5000 is satisfactorily small.
Similarly, the GNSS measurements in 2016 were collected in nighttime by using two multi-constellation Leica GS10 receivers whereas the measurements using the UCM Mekometer ME5000 were carried out during the same days in daytime. No significant discrepancies could also be found among any of the types of distances.
Nevertheless, all these distances constitute but a limited sample so that additional measurements need to be made to draw definitive conclusions and further research has to be conducted in order to assess the validity of this degree of agreement and possibly reduce the standard deviation of the GNSS distances computed.

Conclusions and Future Work
We have proposed and tested new methods to develop a submillimetric GNSS-Based Distance Meter (GBDM) for application to length metrology in open air for distances up to 1 km, as an alternative to more expensive approaches to the problem like two-colour telemeters or laser trackers, which are currently under development. We present the corresponding functional model as well as specific methods to mitigate the main source of error in the GBDM determination of length for baselines up to 1 km: the multipath error. Consistency with the unit of length (SI-metre) has been shown within few tenths of a millimetre by using a calibration baseline of reference having distances traced to the SI-metre and a ME5000 distance meter.
The future application of the GBDM to longer baselines, in the range of 1 to 5 km, can be expected as soon as some difficulties are circumvented: mainly, the correction and corresponding assessment of uncertainty of the impact on the baseline length of the residual tropospheric delay, which should not be neglected especially in the case of relatively large height differences. The evaluation of different absolute antenna calibration models for use in the problem at hand and the corresponding uncertainty assessment is also a question to be addressed in the future.
The present study has been realized with the main focus on the GPS constellation, and the marginal use of satellites from other constellations (Galileo and GLONASS). An optimized, fully multi-constellation, approach needs also to be further developed.