Hybrid Atmospheric Modeling of Refractive Index Gradients in Long-Range TLS-Based Deformation Monitoring

1. Introduction

1.1. Problem Description

Terrestrial laser scanners (TLS) are the active sensors capturing millions of the points per second via evaluating the reflected signal from target surfaces at the TLS. To guarantee the high quality of TLS deliverables (3D point clouds), acquired observations must be assessed in terms of four systematic error sources: instrumental imperfections, atmospheric effects (refraction), target and surface related parameters and scanning geometry (reflectivity). It is reasonable that an understanding of optical effects of the laser, refraction and reflection, is necessary towards designing a more rigorous calibration setup.

Since the emergence of terrestrial laser scanners, they have become crucial in performing engineering geodesy and deformation analyses tasks. For instance, due to safety reasons, the preferred technique of observation for unapproachable areas and high walls, especially in mining and structural sites, is long-range terrestrial laser scanning method, typically within the $500m-1000m$range or more. Long-range scanning demands detailed knowledge of the optical influences on the geometry of the line of sight to genuinely differentiate between the actual geometric displacements and TLS observations uncertainties. Given the optical influences along the line of the sight, geodetic refraction is the dominating effect. Other optical occurrences, especially reflectivity or scattering, are not addressed in this research, but those are the future research topics into the TLS performance.

Geodetic refraction is a deviation of the signal from its direct line of travel, due to the varying velocity of the wave propagating into different mediums of the atmosphere. Therefore, the signal follows the quickest path through the medium to reach its destination (the surface of the targets). A change in refraction might take place numerous times during the period of observations due to variations in atmospheric conditions. The atmospheric variables such as air temperature, atmospheric pressure, relative humidity, etc. are predominantly contributing elements for geodetic refraction [1].

1.2. Significance and Purpose

The fundamental restriction in determination of refraction is non-uniformity of atmosphere conditions (i.e., its turbulence) over the optical path. In addition, no current technological advancement is capable of monitoring the variation of refractive index with respect to the corresponding atmospheric conditions. Previously, observations of atmospheric conditions at both terminals of the sightline were acquired, and the mean calculation of the refractive index at the highest precision of $\left[1-5\right]\times {10}^{-8}$ was employed to correct the measured range, to support millimeter- or sub-millimeter-accuracy of distance measurements over long baselines. For example, range corrections for the refractive index effects $∆r$ over $1000m$will be of the order of sub-millimeter if refractive index is precise at the level of $1\times {10}^{-8}$ ($∆r=1\times {10}^{-8}\times 1000m=0.01mm$). Therefore, ideally, this level of precision or better for refractive index estimation must be satisfied. This principle refers to the second velocity correction of the range measurements [1,2,3].

The current study, for the first time, introduces a generic atmospheric error model for long-range terrestrial laser scanning based on physical refractive index parameterization. Two physical refractive index models - the Ciddor and Closed Formula models - are presented, and further theoretical developments are revised to account for spatial variations of the refractive index - vertically and horizontally. This development results in nonlinear refractive index modelling along the sightline by incorporating varying spatial refractive index gradients within the assumed vertically stratified atmospheric layers, rather than relying on the conventional physical model that applies an average of refractive indices from both endpoints. The robustness of the methodologies is subsequently compared and validated using two geodetic monitoring datasets acquired with two long-range scanners: the Leica ScanStation P50 and the Maptek I-Site 8820. The first dataset, collected at a mine site, enabled long-range scanning exceeding 800 m, while the second, collected at a dam site, allowed extreme vertical viewing angles approaching the zenith (80°).

To achieve high-precision calibration results, both test fields were established within a calibrated network using post-processed GPS control points and onsite terrestrial surveys, delivering $\pm 1mm$ accuracy for the range and $1”$ for the vertical angle observations. Additionally, in-situ atmospheric observation methods were implemented in conjunction with each scan setup. A high precision for refractive index estimation $1\times {10}^{-9}$is sought prior to assessing the accuracy of 3D point coordinates. Accuracy assessments are completed by comparing the residuals between calibrated results and the control points. The results show the improvements in range accuracy - $34\%$ (from $3.6mm$ to $2.4mm$) for the Leica ScanStation P50, and $16\%$ (from $9.2mm$ to $7.7mm$) for the Maptek I-Site 8820 - as well as the improvements in vertical angle accuracy - $44\%$ (from $18”$ to $10”$) for the Leica ScanStation P50 and $20\%$ (from $24”$ to 19”) for the Maptek I-Site 8820. These are accomplished through the implementation of the advanced physical refractive index model, which addresses the non-linearity of refractive index modeling along the line of the sight by incorporating varying gradients of refractive index. Compared to the insignificant improvements obtained by the conventional physical models, these findings are substantial. Consequently, the 3D point coordinates exhibit maximum accuracy improvements of $41\%$ and $18\%$ for the Leica ScanStation P50 and Maptek I-Site 8820, respectively. The associated uncertainties are reduced from the centimeter level to the millimeter level for 3D point coordinates through the implementation of the advanced approach. In addition, it is recognized that points at longer ranges (greater than $\approx 200m)$and/or at steep vertical angles (greater than $\approx 60°)$ benefit more from atmospheric corrections than points located at shorter ranges and/or close to the same horizontal plane as the scanner station.

The hybrid refractive index model, a data-driven approach that integrates neural network techniques with the advanced physical model, ensures further improvements in measurement accuracy - $1.8mm$ and $9″$ for the Leica ScanStation P50, and $2.3mm$ and $15″$ for the Maptek I-Site 8820 - consistently reducing 3D point coordinate errors to a reliable millimeter level. This enhancement is achieved by maximizing the precision of refractive index estimation while minimizing sensitivity to 3D spatial gradients, particularly vertical gradient, thus capturing real-world refractivity patterns along the propagation path. To conclude, the robustness of these findings – based on a comparison of three methods: conventional and advanced physical refractive index models versus the hybrid refractive index model - was further validated across different scanners and varying scanning environments, consistently achieving millimeter-level uncertainty in 3D point coordinates

2. Literature Review

2.1. TLS Principle

TLS is a very high-speed and movable total station which is able to capture millions of points in a second as the consequence of measuring three spherical coordinates, range $r$, vertical angle $v$ and horizontal angle $h$ from a returned signal reflected from a single point. The mathematical conversion from 3D spherical coordinates $\left[\begin{array}{ccc}{r}_p& v_p& h_p\end{array}\right]$into Cartesian coordinates $\left[\begin{array}{ccc}{x}_p& y_p& z_p\end{array}\right]$is applied as below (i.e., $p$ is the number of measured points from $1$ to $n$) [4]:

$${\left[\begin{array}{c}{x}_p\\ y_p\\ z_p\end{array}\right]}_{p=1\dots n}={\left[\begin{array}{c}{r}_p\cos v_p\cos h_p\\ r_p\cos v_p\sin h_p\\ r_p\sin v_p\end{array}\right]}_{p=1\dots n},$$

(1)

Afterwards, the projection is required from scanner coordinate system $j$ to object coordinate system into space $\left[\begin{array}{ccc}{X}_P& Y_P& Z_P\end{array}\right]$using the rotation matrix $M$ and translation parameters of $\left[\begin{array}{ccc}{X}_S& Y_S& Z_S\end{array}\right]$:

$$\left[\begin{array}{c}{x}_p^j\\ y_p^j\\ z_p^j\end{array}\right]=M·\left(\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]-\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]\right),$$

(2)

There are three types of range measurement techniques used by TLS: time of flight (TOF), phase based and waveform digitizer (WFD) (Figure 1). TOF method is the technique that time plays the integral role in capturing the range (the time is recorded between emitting the signal pulse from the instrument and receiving the same pulse at the instrument). TOF uses either pulsed modulation or continuous wave modulation (CW). Whereas phase-based method conveys the data by modulation of the phase of the signal (i.e., the range is determined by the shift in phase between emitted and received signal considering the number of full wavelength). The main distinction between TOF and phase-based principle is that the phase-based technique provides more information to evaluate the entire signal, including signal shape, channel amplification, etc. which enables more accurate distance determination compared to TOF. However, TOF is appropriate for the long-range measurement with the less accurate deliverables due to the expansion of laser spot size within a range. WFD is the advanced combination of both techniques, in which the time between a start and stop pulse is calculated, and the entire received signals is digitized [5].

Additionally, for digital angle recordings of laser scanners, various technologies have been introduced - camera, hybrid and panoramic. Those are classified based on varying horizontal and vertical field-of-view (Figure 2). Then, to determine vertical and horizontal angle for panoramic scanners, following equations must be used to convert the Cartesian coordinates to spherical coordinates: $v=\pi -{\mathit{tan}}^{-1}({\displaystyle \frac{z}{\sqrt{x^2+y^2}}})$ and $h={\mathit{tan}}^{-1}({\displaystyle \frac{y}{x}})-\pi$ [6,7].

In principle, the electromagnetic wave (EM) is the main source of illumination for the TLS measurements (i.e., employed domain for wavelength is typically from visible to far infrared ($400nm$to $2000nm$)). Therefore, the accuracy of observations highly depends on the accuracy of the propagating wavelength and its velocity though the travelling time period. Electromagnetic waves are generally described by wavelength $\lambda$ in $m$, frequency $f$ in $Hz$ or $s^{-1}$ and the propagation velocity $c$ in [1]. The relationship is expressed based on:

$$\lambda ={\displaystyle \frac{c}{f}},$$

(3)

EM in air is influenced by several atmospheric conditions, with the most prominent ones being the variation in air temperature, atmospheric pressure of air, water vapor in air (e.g., humidity), the effects of carbon dioxide content, etc. Those environmental influences change the direction of the line of the sight from the chord (the corrected line of the straight sightline) to the actual ray path, resulting in a range and vertical angle deviation.

Figure 3 shows that the deflection from the chord to actual ray results in different observations on the surface of scanned objects for the range and potentially for the angle observations from the perspective of the apparent ray (i.e., $dv$ is the potential deviated vertical angle). In optics, this phenomenon is defined by Fermat’s principle and Snell’s law as the refraction of a wave. Then, the relationship of the propagating velocity of EM waves into air $c_{air}$compared to the velocity of the identical wave into a vacuum is indicated as the refractive index $n$ (The index is dimensionless).

$$n={\displaystyle \frac{c_{vacuum}}{c_{air}}},$$

(4)

And refractivity $N$ (without the metric unit) can be explained as follows [1,3]:

$$N=(n-1)\times {10}^6,$$

(5)

The type of refraction shown in Figure 3 is a convex condition. When the refractive index from the first to the second medium drops ${(n}_2<n_1)$, it refers to a convex condition, otherwise if the refractive index rises from the first to the second medium $(n_2>n_1)$, the refraction condition is concave [10].

2.2. Review and Results on Mathematical Developments of Refraction

The difficulty of mitigating the relevant impact of refractive index is non-uniformity of atmospheric conditions (i.e., the turbulence of refractivity) over the optical path. To deal with this limitation, two methodologies were presented in the literature: direct method, or indirect method. The direct method of refractive index observations is implemented using interference refractometer [11]. This technique brings some disadvantages despite the direct measurement of refractive index. One of the disadvantages is that the refractive index is determined independent of atmospheric variations of the air (i.e., it does not physically reflect the real condition of atmosphere). Moreover, the internal calibration is required to align the initial crude measurement relative to the length of the refractometer cell, that introduces more computational efforts [11]. The alternative method is the indirect method. In this technique, the refractive index is derived from measurement of atmospheric parameters. Thus, to achieve the optimum performance of the index, it is recommended that all atmospheric conditions are precisely obtained at both terminals of the line, and the mean calculation for refractive index is determined at least at the precision level of $\left[1-5\right]\times {10}^{-8}$to guarantee millimetre- or sub-millimetre-precision range [1].

Mathematical developments of refractive index with respect to the atmospheric variables have been updated regularly since the 19th century. The index has been repeatedly studied and investigated by different scientists to achieve better results (i.e., the studies are listed by years [12]; [13,14]; [15]; [16]; [17,18]; [19]; [20,21]; [22,23]; [11], the model adopted by International Association of Geodesy called Closed Formula in 1999 [24]; [25], and [26]). Under each notion, at least three elements of atmospheric conditions have been considered - air temperature, atmospheric pressure and relative humidity.

M. Sabzali et al [27] conducted extensive research on the proposed models of refractivity and investigated the impacts of refraction on the simulated and real datasets for the range observations. The aim of the work was to update atmospheric modelling for the range measurements. Their findings were validated over existing range measurement techniques: TOF and phase shift principles. After improving the atmospheric error model, the variation of reference index $dn$ over the entire line of path is verified. In Equation 6, $dn$ is the comparison between two actual atmospheric conditions (obtained from two varying mediums):

$$dn\times {10}^6=-0.93dT+0.28dP-0.039de,$$

(6)

Where, each change in air temperature either $dT\left(K\right)$ or $dt(℃),$ atmospheric pressure $dP$ ($hPa$) and partial water vapour pressure (humidity) $de$ ($hPa$) respectively introduces $-0.93ppm,+0.28ppm$ and $-0.039ppm$ refraction over the measured range $r_o$ [24].

The impact of refractivity on the angle measurements - vertical and horizonal angle - were also separately studied (Figure 4).

The improved error models are modified for refracted vertical angle $dv$ and refracted horizontal angle $dh$ based on the experimental analysis in [29,30]:

$$dv=\left[0.00026dT-0.00004dP\right]·r_o,$$

(7)

$$dh=\left[0.001\times {10}^{-6}dT+1.7\times {10}^{-6}dP\right]·r_o,$$

(8)

Here, one single unit rise in temperature and pressure between the mediums respectively leads to $0.26”$ and $-0.04”$ refraction in measured vertical angles within the assumed range observation of $1000m$. However, the variations for the horizontal angles are considerably smaller (i.e., given insignificant effect of humidity on angle measurements). Further, M. Sabzali. et al [30] acknowledged that the effects of refraction on vertical angles for points located close to nadir and zenith is maximum due to the sensitivity between vertical and horizontal gradients of refractions.

3. Methods

The quantification of refractive index $n$ along the line of the sight must be precisely addressed when the beam ray experiences different refraction conditions with respect to the atmospheric variations in each medium. This might occur several times during time of the measurement. Therefore, the error model for optimization of the systematic error caused by atmospheric effects, from Equation 1, can be rewritten as follows:

$${\left[\begin{array}{c}{x}_p\\ y_p\\ z_p\end{array}\right]}_{p=1\dots n}={\left[\begin{array}{c}({r+dr)}_p\cos ({v+dv)}_p\cos {(h+dh)}_p\\ {(r+dr)}_p\cos ({v+dv)}_p\sin ({h+dh)}_p\\ {(r+\mathit{dr})}_p\sin ({v+dv)}_p\end{array}\right]}_{p=1\dots n},$$

(9)

where, $\left[\begin{array}{ccc}dr& dv& dh\end{array}\right]$ refers to the refracted range, refracted vertical angle and refracted horizontal angle, respectively. To reduce the complexity of the problem, it was suggested that accurate observations of atmospheric conditions at both terminals of the sight line are measured, and then the mean calculation of the refractive index is applied as the real refractive index [1]. These correction factors are referred to as the first velocity correction, and its comparison with the reference refractive index $n_{ref}$ is the second velocity correction for range measurements. The second velocity correction is typically implemented in order to ensure the required precision for estimation of refractive index (at least better than $[1-8]\times {10}^{-8}$). Note, in most cases, for either end, the defaulted value within the instrument is assumed as the reference refractive index [3,31]. Since TOF is the scanning mechanism behind long-range scanners for range measurements, the first and second velocity correction factors must be simultaneously substantiated after precise observations of atmospheric variables, through a physical refractive index model.

As discussed earlier, there are a number of physical refractive index models. Among all sets of calculations, Ciddor’s parameterisations provide more robust results than the previous versions [22,23]. Under this formulation, it was noticed that this setup is appropriate over a broader range of wavelengths (within $300nm-1690nm$). Additionally, Ciddor physical model supports more flexibilities under extreme environmental conditions [26,27]. To implement Ciddor’s refraction model, the three following steps must be implemented:

1.

The first step is to differentiate the phase refractive indices1 of standard air2 $n_{st}$ and water vapor $n_{wv}$ as the function of the wavelength and irrespective of atmospheric variables:

$${N_{st}=\left(n_{st}-1\right)\times 10}^8={\displaystyle \frac{a_1}{a_2-{\sigma }^2}}+{\displaystyle \frac{a_3}{a_4-{\sigma }^2}},$$

(11)

Here, the wave number $\sigma$ in ${\mu m}^{-1}$is reciprocal of wavelength $\lambda$ ($\mu m)$.

The empirical coefficients are listed in Table A1 in the Appendix A.

$$N_{wv}={\left(n_{wv}-1\right)\times 10}^8=cf\left(a_5+a_6{\sigma }^2-a_7{\sigma }^4+a_8{\sigma }^6\right),$$

(12)

The correction factor $cf$= $1.022$ is considered (dimensionless).

The group refractive indices are also computed as below:

$${N_{st}=\left(n_{st}-1\right)\times 10}^8={\displaystyle \frac{a_1\left(a_2+{\sigma }^2\right)}{{\left(a_2-{\sigma }^2\right)}^2}}+{\displaystyle \frac{a_3(a_4+{\sigma }^2)}{{(a_4-{\sigma }^2)}^2}},$$

(13)

$$N_{wv}={\left(n_{wv}-1\right)\times 10}^8=cf\left(a_5+a_9{\sigma }^2-a_{10}{\sigma }^4+a_{11}{\sigma }^6\right),$$

(14)

2.

Next step is to compute the refractive index based on the atmospheric conditions, density components of the dry air ${\rho }_{axs}$ and ${\rho }_a$, and the moist air ${\rho }_{wv}$and ${\rho }_w$ with corresponding values of compressibility of air COM:

$${\rho }_{axs}={\rho }_{wv}={\displaystyle \frac{P{M}_a}{\left(COM\right)RT}}\left(1-x_w\left(1-{\displaystyle \frac{M_w}{M_a}}\right)\right),$$

(15)

where,$T,$ is temperature $K,$ $P$ is the pressure $\left(hPa\right)$, $x_w$ is the water vapor pressure component of the air depending on the humidity $\left(hPa\right)$ (as three major atmospheric components for this contribution), $M_a$is the molar mass of water vapour containing $x_c$ $ppm$ of ${CO}_2$ $kg/mol$, $M_w$is the molar mass of water vapor $(=0.018015kg/mol)$, and $R$ is the gas constant ($=8.314651$ $J{mol}^{-1}{K}^{-1})$. Then, compressibility $COM$ based on each air conditions - either standard dry air or pure water vapor - are computed:

$$COM=1-\left({\displaystyle \frac{P}{T}}\right)\left[a_{12}+a_{13}t+a_{14}{t}^2+\left(a_{15}+a_{16}t\right)x_w+\left(a_{17}+a_{18}t\right)x_w^2+{\left({\displaystyle \frac{P}{T}}\right)}^2\left(a_{19}+a_{20}{x}_w\right)\right],$$

(16)

$${\rho }_w={\displaystyle \frac{P{M}_w{x}_w}{\left(COM\right)RT}},$$

(17)

$${\rho }_a={\displaystyle \frac{P{M}_a{(1-x}_w)}{\left(COM\right)RT}},$$

(18)

Here, $t$ is the temperature in $℃\left(t=T-273.15\right).$

3.

Ultimately, the combined evaluation of both refractive indices - under dry air and water vapor component - is determined by [22,23]:

$$N_{ciddor}=(n_{ciddor}-1){\times 10}^8=\left({\displaystyle \frac{{\rho }_a}{{\rho }_{axs}}}\right)\left(n_{st}-1\right)+\left({\displaystyle \frac{{\rho }_w}{{\rho }_{wv}}}\right)\left(n_{wv}-1\right),$$

(19)

Ciddor’s parameterisations have been later adopted by the International Association of Geodesy (IAG) in 1999, as the standard equation for calculating index of refraction for geodetic instruments operating within the visible and near infrared waves [24]. The principle is referred to as the Closed Formula model. IAG’s proposal provides more accurate results under more extreme temperature, pressure and humidity conditions through simplification of the Ciddor’s principles and less computational skills. Thus, the methods to achieve the group refractive index $n_g$as a function of $\lambda$ in $\mu m$ is straightforward as follows:

$$N_g=\left(n_g-1\right)\times {10}^6=b_1+{\displaystyle \frac{b_2}{{\lambda }^2}}+{\displaystyle \frac{b_3}{{\lambda }^4}},$$

(20)

Afterwards, the group refractive index $n_{iag}$under either standard air condition or water vapour is computed:

$$N_{iag}=(n_{iag}-1)\times {10}^6={\displaystyle \frac{b_4{N}_{g}P}{b_{5}T}}-{\displaystyle \frac{b_{6}e}{T}},$$

(21)

All empirical coefficients are listed in Table A1 in the Appendix A.

Given each model, either the Ciddor or Closed Formula, three spatial variations of the refractive index $n\left(x,y,z\right)$in 3D Cartesian coordinates can be parameterised as the gradient of the refractive index $\nabla n$:

$$\nabla n={\displaystyle \frac{\partial n}{\partial x}}\overrightarrow{i}+{\displaystyle \frac{\partial n}{\partial y}}\overrightarrow{j}+{\displaystyle \frac{\partial n}{\partial z}}\overrightarrow{k},$$

(22)

The elements on the horizontal plane $\left[\begin{array}{cc}\partial n/\partial x& \partial n/\partial y\end{array}\right]$ affecting the horizontal directions are called horizontal gradients of refraction, while $\partial n/\partial z$refers to the vertical gradient of refractive index impacting the vertical directions. Therefore, the gradient of refractive index is rewritten as the function of the atmospheric variables:

$$\nabla n=\left({\displaystyle \frac{\partial n}{\partial T}}\nabla T+{\displaystyle \frac{\partial n}{\partial P}}\nabla P+{\displaystyle \frac{\partial n}{\partial e}}\nabla e\right),$$

(23)

where, $\nabla T=\left[\begin{array}{ccc}{\displaystyle \frac{\partial T}{\partial x}}& {\displaystyle \frac{\partial T}{\partial y}}& {\displaystyle \frac{\partial T}{\partial z}}\end{array}\right](K/m$ or $\mathrm{℃}/m),\nabla P=\left[\begin{array}{ccc}{\displaystyle \frac{\partial P}{\partial x}}& {\displaystyle \frac{\partial P}{\partial y}}& {\displaystyle \frac{\partial P}{\partial z}}\end{array}\right](hPa/m)$ and $\nabla e=\left[\begin{array}{ccc}{\displaystyle \frac{\partial e}{\partial x}}& {\displaystyle \frac{\partial e}{\partial y}}& {\displaystyle \frac{\partial e}{\partial z}}\end{array}\right](hPa/m)$ are horizontal and vertical gradient of air temperature, atmospheric pressure and humidity of air, respectively [32].

To investigate different gradient components of the refractive index, the stable stratification of the atmosphere - based on air temperature - is sometimes supposed in analytical studies. The stable stratification condition is defined when air temperature decreases gradually with height. In contrast, unstable stratification occurs when the air temperature decreases rapidly with height, or changes irregularly with height. The vertical temperature gradient is generally described as the variation of temperature vertically under stable stratification conditions. In the following studies [29,33], different layers of the atmosphere - with respect to the height from ground level, from $0m$ (directly above the ground) to over $100m$ - were categorised, and the corresponding practices for computing the vertical temperature gradient were presented. In summary, vertical temperature gradients are substantially intense in the layers close to the ground surface within a range of $0$ to $3m$ (between $-47K/m$ and $+20K/m$) and drops to the small value $-0.006K/m$ at the highest assumed level.

Alternatively, the variation of temperature horizontally or laterally, impacting the horizontal direction of refractive index is the horizontal temperature gradient. There have been several methodologies to determine the horizontal temperature gradient. For example, B. G. Bomford [2] assumed that for an approximate distances of $1000m$ horizontally and $3000m$ vertically, with $5K$ temperature rise between the terminals, the difference in average horizontal temperature gradient is nearly $0.005K/m$ (i.e., the increase of $5K$ per kilometer in lateral difference). Another theory in the lowest layer of atmosphere – where a line of sight is two meters above the ground level with the vertical temperature gradient of $0.3K/m$ - the horizontal temperature gradient is checked, and it is assumed negligible.

Comparing these two gradients, the effect of atmospheric variables horizontally is trivial. Therefore, Equation 24 depicts the relationship between vertical gradient of refractive index ${\displaystyle \frac{\partial n}{\partial z}}$and the vertical temperature gradient ${\displaystyle \frac{\partial T}{\partial z}}$ as well as the vertical pressure gradient ${\displaystyle \frac{\partial P}{\partial z}}$and vertical humidity gradient ${\displaystyle \frac{\partial e}{\partial z}}$ [27]:

$${\displaystyle \frac{\partial n}{\partial z}}=\left[{\displaystyle \frac{11.27e}{T}}-{\displaystyle \frac{76P}{T^2}}\right]{\displaystyle \frac{\partial T}{\partial z}}+{\displaystyle \frac{76}{T}}{\displaystyle \frac{\partial P}{\partial z}}+-{\displaystyle \frac{11.27}{T}}{\displaystyle \frac{\partial e}{\partial z}}$$

(24)

As discussed earlier, the primary contributor to the vertical gradient of refraction is the vertical temperature gradient, with minimal influence from other vertical gradients. Further developments regarding the pressure and humidity vertical gradients are presented in Appendix B.

Consequently, under these assumptions, the corrected range $r_c$ and corrected vertical angle $v_c$ can be expressed in terms of the refracted range $dr$ and refracted vertical angle $dv$, respectively. These are derived by integrating of the refractive index effects over the entire length of actual ray (observed range $r_o$) (Figure 1):

$$r_c={\int }_o^{r_o}{\displaystyle \frac{dr}{n\left(z\right)}},$$

(25)

$$v_c=dv+v_o={\int }_o^{r_o}{\displaystyle \frac{1}{n}}{\displaystyle \frac{dn}{dz}}dr+v_o,$$

(26)

Also, a refracted horizontal angle $dh$ is reparametrized in terms of horizontal gradient of refractive index ${\displaystyle \frac{dn}{dy}}:$

$$dh={\int }_o^{r_o}{\displaystyle \frac{1}{n}}{\displaystyle \frac{dn}{dy}}dr,$$

(27)

Accordingly, to achieve an optimal performance of the refractive index correction through a physical model, the refractive index is approximated by the averaged values obtained at both terminals of the line of sight $z_0$and $z_1$ [1]. This means that refractive index is expressed as (Section 5.1):

$$\mathit{n}\left(\mathit{z}\right)\approx {\displaystyle \frac{1}{2}}\left[\mathit{n}\left({\mathit{z}}_0\right)+\mathit{n}\left({\mathit{z}}_1\right)\right],$$

(28)

However, the imposed simplification underestimates the potential non-linear vertical variations of the atmosphere along the actual ray path, as the laser beam experiences multiple refraction abnormalities when traversing different atmospheric layers. Therefore, in this research, a more accurate physical model is proposed to account for the varying vertical gradients of the refractive index along the propagation path – referred to as advanced physical refractive index model. According to Equation 25, it can be rewritten in the following expression:

$$r_c={\int }_o^{r_o}{\displaystyle \frac{dr}{n\left(z\left(r\right)\right)}},$$

(29)

where $z\left(r\right)$represents the height profile along the laser path, with$r$ being the counter that moves between $0$ and $r_o$. The refractive index is approximated as the nonlinear function to the range rather than applying the mean value between two refractive indices:

$$\mathit{n}\left(\mathit{z}\left(\mathit{r}\right)\right)\approx \mathit{n}\left[{\mathit{z}}_0+{\displaystyle \frac{{\mathit{z}}_1-{\mathit{z}}_0}{{\mathit{r}}_{\mathit{o}}}}\mathit{r}\right],$$

(30)

Therefore, from Equation 26, the refracted vertical angle $dv$ can be (Section 5.2):

$$\mathit{d}\mathit{v}={\int }_{\mathit{o}}^{{\mathit{r}}_{\mathit{o}}}{\displaystyle \frac{1}{\mathit{n}\left(\mathit{z}\left(\mathit{r}\right)\right)}}{\displaystyle \frac{\mathit{d}\mathit{n}}{\mathit{d}\mathit{z}}}\mathit{d}\mathit{r},$$

(31)

Note, in either case, to guarantee millimetre- or sub-millimetre-accuracy for the observations, a precision of at least $\left[1-5\right]\times {10}^{-8}$must be accomplished for the estimation of refractive index [1]. Then, it enables reducing sensitivity to the spatial variations in the refractive index and enhancing the overall accuracy of 3D point coordinates.

To represent refractivity along the entire line of sight, the illustrated techniques on physical model establish a reasonable relationship between environmental parameters (e.g., in-situ atmospheric recordings) and wave number (i.e., scanner wavelength) as the inputs, and the real-world refractivity along the beam path as the output. However, to achieve rigorous precision and consistency in TLS-based physical refractive index modeling, a hybrid physical-data-driven model is proposed as a follow-up. This hybrid model integrates results from the advanced physical model with a neural network approach, and its function is justified through field-validated outcomes from previously established physical models. It guarantees optimal millimeter-level accuracy in 3D point coordinates and enables consistency checks across two long-range scanners under varying atmospheric conditions (Section 5.3).

In short, a neural network is a machine learning algorithm simulated from the structure of the human brain. It generally consists of a variety of layers from interconnected neurons, which each one collects the inputs and interacts with the result to the next consecutive layers to generate rigorous outputs. These neurons utilize certain mathematical algorithms and adjust internal weights during each training interval to minimize the prediction errors for both datasets, according to the received residuals [34]. The minimization of residuals of 3D spherical coordinates signifies the lowest ultimate uncertainty for the 3D point coordinates (on the order of millimetre or sub-millimetre relative precision). It ultimately ensures the robust prediction by the most accurate output - refractive index and its spatial gradients along the laser path.

In summary, three different methods are tested to improve calibration uncertainty (Table 1). To achieve high-precision calibration setups, two geodetic test fields, mine site and dam site, were established within a calibrated network using post-processed GPS control points and onsite terrestrial surveys, delivering $\pm 1mm$ accuracy in range and $1”$ in vertical angle observations. In addition, in-situ atmospheric observations were collected at each scan station to improve refractive index modeling. Using the proposed approaches, refractive index estimation with high precision is implemented prior to 3D point coordinate accuracy assessments.

4. Data Experiments

The above theoretical developments were tested on real case studies acquired from a mine site and a dam site (Figure 5). The mine site experimental test field examines long-range scanning with the maximum range of $846.304m$, while the dam site experimental test field provides the flexibility for investigation of a steep vertical angle from the bottom of the dam to the dam crest (maximum vertical angle captured on-the-site is $80°4^{’}22”$). At the mine site, the data field capture was set up within a calibrated network using eight GPS control points distributed at different elevations across the site (red points shown in Figure 5). The reason for distributing the control points at varying heights is to investigate the varying vertical gradient of refractive indices across different horizontal stratifications of the atmosphere for the advanced hybrid model (e.g., $74.936m$for Station 1, $84.803m$for Station 2, and $128.531m$for Station 3 (Table A2 in Appendix C)). At the dam site, a Leica Nova MS60 MultiStation was used to measure $14$ black and white targets established on the semi-vertical dam walls (i.e., range and angular accuracies of the Leica MS60 are $1mm+1.5ppm$ and $1”$, respectively3).

Furthermore, GPS control points for the mine site were collected on-the-site using static mode and have been post processed after the field collection at the office to reduce to the highest accuracy within $1$ to $5mm$. GPS control coordinates and survey control marks for both datasets are listed in Table A2 and Table A3 in the Appendix C, respectively. For scanning, two long-range scanners - Leica ScanStation P50 (measurement range $1000m$) and Maptek I-Site 8820 (measurement range $2000m$) were employed (Figure 6), and Table 2 indicates technical specifications of the scanners, reported by the manufacturers, and contains the scanning characteristics used in this research.

The datasets from three nominal scanner stations were captured under identical filed instructions on the 10th of December 2024 (mine site) and the 15th of February 2025 (dam site), during working hours from $8:00$ to $17:00$. During scanning, long-range mode within the scanners was activated, and all default correction factors including instrumental atmospheric refraction were switched off. At least two scans from each station with the maximum possible instrumental resolution were acquired (Table 2) (i.e., each with different horizontal orientation).

The environmental conditions of the sites were precisely recorded during scanning time using a Kerstal 2500 weather meter sensor. The reported precisions for the temperature and pressure are $0.5℃$ and $1.5hPa$, respectively. The link regarding technical provisions of the thermometer was provided6. Figure 7 compares two different temperature recordings across two geodetic sites.

To optimally comprehend the level of atmospheric variations across the field, the atmospheric recordings were achieved by at least two weather meter sensors, and the measurements were initiated from the surrounding of each scanner station and continued onsite close to every target location – from the nearest to furthest target location - during scanning time. For example, at $9am$ when the scanner was set up at the first station, the observed temperature close to the station was $43℃$ at the mine site. Across the entire test field, the temperature varied by $\pm 2°C$ relative to the station’s record (i.e., the atmospheric data attached to each scanner station). By $3pm$, however, the temperature had increased to $46°C$ across the whole site. In addition to temperature recordings, the atmospheric pressure of both sites was observed using the same sensor. Compared to temperature, the atmospheric pressure across the areas was considerably stable during the observation time ($1009hPa$and $1012hPa$ for the mine and dam site, respectively).

5. Data Analysis and Discussions

The datasets from both scanners were collected in the calibrated test field under uncontrolled environmental conditions. For the mine site dataset, the range consistency calibration method is implemented, while for the dam site dataset, on-the-site calibrated angle measurements using terrestrial surveying are undertaken. In short, the range consistency method provides the geometric accuracy as the result of verifying the ranges between each pair of corresponding control points remain invariant across multiple scanner stations [35,36,37]. The existence of eight control points for the mine site dataset delivers $28$ Euclidean ranges for each station, and 14 targets for the dam site dataset offer the same number of angles per scanner setup. Total redundancy for each dataset is six times those numbers (depending on the number of scans). Then, the following step is the accuracy assessment of different proposed physical models in long-range terrestrial laser scanning.

After importing scanned data for each instrument separately using Maptek PointStudio 2024.1.17 and Leica Register 360 software8, each of the corresponding control points was manually selected, and spherical coordinates using the developed MATLAB codes were determined. There was no software registration implemented for the calibration arrangement. The reason is that either software registration (automatic point-to-point or cloud-to-cloud) or manual registration imposes an additional root mean square error (RMSE), originated from the software comparison between the clouds. Importantly, the in-situ atmospheric recordings were attached to each scanner station (Figure 7).

Primarily, the reliability of the datasets is examined using hypothesis tests. The objective is to eliminate the measurements containing outliers due to the manual target selection and reducing the potential propagation of noise, which might otherwise degrade the precision of refractive index estimation, leading to lower accuracy of 3D point coordinates. The hypothesis test compares the weighted $W$sum of the squares of the residuals $V^{t}WV$ against the chi-square distribution${\chi }_{\alpha ,r}^2$, with redundancy numbers $r$, and the significance level $\alpha$. Given the assumed significance level at least $2\%$ for both datasets and corresponding redundancy, the test fails if $V^{t}WV$ is greater than critical value of the distribution (outliers exist in the measurements), or if this is smaller than the critical value, the test passes (better precision than prescribed) [38] (where, $v$ and $\sigma$ are residuals and standard deviation, respectively).

$$V^{t}WV={\displaystyle \frac{v_i^2}{{\sigma }_i^2}}+\dots .+{\displaystyle \frac{v_n^2}{{\sigma }_n^2}},$$

(32)

The underlying assumption of the hypothesis test is that outliers can be detected in a reasonable manner (internal reliability), and the impacts of other undetected outliers are insignificant (external reliability). Then, given either physical refractive index model (the Ciddor or the Closed Formula model), the refractive index is determined at the maximum level of precision, optimally reflecting the real-world atmospheric conditions along the sightline. Subsequently, using Equations either 25, 26 or 28, the estimated precision of the refractive index directly affects the accuracy of the 3D spherical coordinates. The uncertainty of the 3D Cartesian coordinates $\left[\begin{array}{ccc}{\sigma }_x& {\sigma }_y& {\sigma }_z\end{array}\right]$ is evaluated using the principle of propagation error, representing the posteriori accuracy of the advanced model - whether physically or hybrid:

$$\left[\begin{array}{c}{{\sigma }_x}^2\\ {{\sigma }_y}^2\\ {{\sigma }_z}^2\end{array}\right]=\left[\begin{array}{c}{(\cos v\cos h{\sigma }_r)}^2+{\left(-\mathrm{r}sinv\cos h{\sigma }_v\right)}^2+{(-\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}v\sin h{\sigma }_h)}^2\\ {\left(\cos v\sin h{\sigma }_r\right)}^2+{\left(-r\sin v\sin h{\sigma }_v\right)}^2+{\left(r\cos v\cos h{\sigma }_h\right)}^2\\ {(\sin v{\sigma }_r)}^2+{(r\cos v{\sigma }_v)}^2\end{array}\right],$$

(33)

Here, $\left[\begin{array}{ccc}{\sigma }_r& {\sigma }_v& {\sigma }_h\end{array}\right]$ represent the posteriori accuracy of the 3D spherical coordinates. However, the accuracy assessment is interpreted as relative precision, given the accuracy of the control network $(\pm 1\sigma =1mm$ and $1”$). Note that identical data analysis procedures are followed for both datasets. Figure 8 highlights a broad summary of the proposed calibration methodologies, considering the mentioned criteria.

5.1. Physical Refractive Index Model: Conventional Approach

The mine site dataset examines long-range scanning within a calibrated network through a distribution of eight control points. To initiate the data analysis, the computed inter-target ranges from selected targets in each station’s point cloud are determined and validated against control points through a range consistency method (i.e., a total number of 168 ranges for the entire network). The range consistency calibration method ensures that the ranges between each pair of corresponding control points in the scanned data remain consistent across different scanner stations [35,36,37].

For pre-processing of the physical model, reliability test on selected targets is recognized (using the hypothesis test Equation 32), resulting in 18 measurements for Leica ScanStation P50 and 12 measurements for Maptek I-Site 8820 being detected as the outliers and eliminated from the dataset. Figure 9 shows the results of comparing each pair of control ranges with the average residuals, considering their height differences (i.e., the average ranges are derived from observations made at multiple scanner stations).

The comparison above underlines that after applying the lowest significance level for hypothesis test ($2\%$), the residuals experience the noticeable variations. However, the obtained standard deviations of the measured ranges are $3.6mm$ for the Leica ScanStation P50 and $9.2mm$ for the Maptek I-Site 8820. Note, with the Maptek device, the larger standard deviation had been anticipated since the scanner has been operated in open pit mining sites and has not recently been under the uncertainty assessment. Additionally, the nonlinear relationship between the ranges and the residuals, for both scanners, highlights the complexity of atmospheric modelling for long-range terrestrial laser scanning. To date, conventional refractive index correction models for range measurements have typically been assumed to be linear, based on the relative consistency of the refractive index along a horizontal path, using the mean value from both ends of the sightline (Equations 25 and 28) [1,27].

The following reasons elaborate on the limitations of the conventional approach of the physical model. Firstly, vertical stratifications of atmospheric air conditions are substantiated and play a vital role in the atmospheric correction modelling of terrestrial laser scanners [29,33]. M. Sabzali et al. [29] introduced different horizontal layers for the atmosphere depending on the height from ground level, according to the vertically stratified atmosphere. These layers commence from the layer closest to the ground surface, from zero to $3m$ (lowest layer, with the most influential variations in terms of atmospheric conditions), $3-20m$ (intermediate layer), $20-100m$, and above $100m$ (highest layer, with the least influence in terms of atmospheric variations) (Figure 6). It is acknowledged that the vertical variation of the atmosphere does not remain stable for the laser path passing within those layers, particularly for vertical temperature gradients. This vertical stratification of the atmosphere is a key consideration for improving current TLS atmospheric correction models. Due to the extended vertical field of view of terrestrial laser scanners (Table 1), a large number of points are observed close to or at the zenith and nadir when the vertical field of view is maximized, dissimilar to the terrestrial surveying in which more points are captured through the limited vertical field of view near the horizon. Therefore, the observed ranges from terrestrial laser scanners experience varying atmospheric gradients along the path - the varying vertical gradient of refractive index - impacting not only the range but also the vertical angle at different horizontal layers of the atmosphere.

Figure 6. Different conditions of refractive index modelling in TLS. Note, scale is not preserved in this figure (The list of control points is available in Table A2 and Table A4 in Appendix C).

Figure 6. Different conditions of refractive index modelling in TLS. Note, scale is not preserved in this figure (The list of control points is available in Table A2 and Table A4 in Appendix C).

Thus, points at different height levels receive variable residuals, regardless of their measured ranges. For instance, three peak residuals correspond to different height differences and range observations (i.e., $-2.846m$for the $r_{18}=153.915m$, $9.866m$for $r_{12}=432.298m,$and $53.595m$ for $r_{13}=425.325m$). The first two are located within short- and mid-range observation, but their height differences indicate that both terminals of these lengths lie within the lowest atmospheric layer ($0$ to $3m$), and between the lowest and the intermediate atmospheric layers, respectively ($3-20m$). Whereas, for the latter pair of points, the travelling path extends from the lowest to the highest atmospheric layer. Hence, the sightline travels from the layers close to the ground, with substantially noticeable atmospheric variations, to the highest layer with considerably insignificant atmospheric influences (Figure 6, condition (a)).

In contrast, each pair of points at the highest layer or above (Figure 6, condition (b) and (c)), regardless of their elevation and range, is less affected by atmospheric components, such as the ranges $200.091m,354.867m$, $577.662m$, $636.217m$, $652.411m$, and $846.304m$ encompassing varying height differences (Table A4 in the Appendix C). The minimum, maximum, and negligible cases, as shown in Figure 6, represent the bounds of refractive index correction that must be applied in the advanced model.

Given the arguments above, Z-coordinates of points observed at greater vertical viewing angles relative to the station require more significant atmospheric corrections than those of the X- and Y-coordinates at similar angles (e.g., Figure 6, the maximum cases). Nevertheless, for the points observed at smaller vertical angles - such as those lying close to same horizontal atmospheric layer as the scanner station - the atmospheric correction is often negligible (e.g., Figure 6, negligible cases). This simplification leads to an underestimation of atmospheric effects in conventional surveying tasks.

To investigate the introduced sensitivity in the vertical direction, the second dataset from the dam site provides significant vertical angle variations from three nominated scanner stations - from the bottom of the dam wall (lowest layer of atmosphere) to the survey targets on the dam wall and the dam crest (highest layer of atmosphere). In this experiment, 14 survey targets were positioned within a shorter range of scanning but at the higher vertical viewing angle (ranging from $45°$ to $80°$). The spherical coordinates were determined and validated against on-site survey data (i.e., in total, 84 vertical angles were derived for the entire network). After identifying seven and eleven outliers for the Leica ScanStation P50 and Maptek I-Site 8820, respectively, the results are depicted to compare the on-site surveyed vertical angles with the averaged Z-coordinate residuals (extracted from the average of measured coordinates), considering their measured ranges (Figure 7).

Figure 7. Control vertical angles, the ranges and the Z residuals (Dataset 2: dam site).

Figure 7. Control vertical angles, the ranges and the Z residuals (Dataset 2: dam site).

The large vertical angles are associated with increased Z residuals (i.e., the standard deviations of the measured vertical angles are $18”$ for the Leica ScanStation P50 and $24”$ for the Maptek I-Site 8820). This outlines that observations at vertical angles beyond $\pm 50°$ exhibit greater sensitivity to atmospheric refraction effects than those within $\pm 50°$ (such as, Dataset 1: mine site). Figure 6 illustrates the maximum cases of vertical gradient of refractive indices under condition (a), highlighting variations from the lowest to either the intermediate or the highest atmospheric layers.

In summary, when both ends of the sightline are located within the intermediate or higher atmospheric layers, refraction effects remain stable in relatively same horizontal layer (negligible cases; Figure 6, conditions (b) and (c)). However, when either end of the sightline is placed within the intermediate or the lower atmospheric layers, the impact of atmospheric refraction becomes more significant and variant with respect to the vertical gradient of refractive index (maximum cases; Figure 6, conditions (a) and (c)). Since terrestrial laser scanners typically transmit signals from one end of the sightline, the systematic error in relation to refractivity increases with the rise in the vertical angle from the horizontal plane in a nonlinear manner. This comparative discussion underscores the requirement for a physical model to account for the variability of temperature gradient along the laser path. The advanced physical model addresses this by incorporating varying vertical temperature gradients (Equations 29 and 31) for the points observed with a long-range baseline (greater than $200m$) and a steep vertical angle (larger than $\pm 50°$) [27]. Consequently, the interdependence of spatial gradients of refractive index can be resolved through the high-precision determination of the refractive index along the travelling path [28].

5.2. Physical Refractive Index Model: Advanced Approach

To overcome the limitations of the conventional physical model, the advanced physical model was introduced. The conventional approach of physical model, which assumes a uniform refractive index based on the average of two terminals, underestimates refraction impacts along the path with varying vertical gradients of refractive index, particularly at long ranges and steep vertical angles. The advanced model aims to incorporate vertical variations in refractive index, thus accounting for the nonlinear influence of refractivity along the propagation path (Equations 29 and 31). In practice, the model assumes different segmentations of a straight-line sight path through different atmospheric layers, with each layer characterized by its corresponding vertical temperature gradient and vertical refractive index gradient. This approach depicts the vertical stratification of the atmosphere more realistically, especially for the TLS observations near the zenith and nadir where sensitivity to vertical temperature gradients is maximal. Accordingly, the required vertical temperature gradients for each layer are summarized in Table 3. Using Equation 24, the corresponding vertical gradient of refractive index can be computed. These provide the foundation for improved TLS atmospheric correction – referred to as the advanced refractive index model.

By integrating these vertical gradients along the propagation path, the advanced model provides corrected ranges and vertical angles that more accurately reflect real-world atmospheric conditions, bridging the gaps of the conventional approach. The results from the advanced model are then compared with those obtained from the conventional physical model (i.e., priori residuals represent the computed residuals after applying the conventional Ciddor refractive index model, while posteriori residuals correspond to the computed residuals after implementing the advanced Ciddor refractive index model (Figure 8)). Table 4 also compares the priori and posteriori uncertainties of the range observations in the range consistency method, expressed as root mean square error (RMSE) values. These highlight the improvements accomplished through the advanced model over the physical model.

Figure 8. Comparison of range residuals in the range consistency method for the Leica ScanStation P50 ($m)$ (Dataset 1: mine site)9.

Figure 8. Comparison of range residuals in the range consistency method for the Leica ScanStation P50 ($m)$ (Dataset 1: mine site)9.

Results indicate that the conventional physical refractive index model, based on average refractive indices at both endpoints of the range, results in consistency of the refractive index correction along the path, mostly account for the horizontal layer of the atmosphere, and is not able to significantly reduce range accuracy over the long baselines. As discussed earlier, this consistency remains important and cannot be overlooked due to the sensitivity of range observations to 3D spatial gradients of the refractive index (impacting the 3D spherical coordinates). However, a more advanced refractive index model, which accounts for the height profile within each vertically stratified atmospheric layer, less sensitive to the spatial gradients, leads to improved accuracy by approximately $34\%$ and 16% for the Leica ScanStation P50 and Maptek I-Site 8820, respectively. Note, both physical models provide identical RMSE values (priori uncertainties). This confirms the suitability of both physical models for TLS atmospheric correction.

Moreover, the moderate improvements in the range consistency method does not necessarily guarantee an identical level of improvements in 3D point coordinates accuracy. Due to the existence of large baselines and limited vertical viewing angles (i.e., the maximum vertical angle in mine site is approximately 14°), the sensitivity between vertical gradient of refractive index ${\displaystyle \frac{dn}{dz}}$and range refractive index $n$ cannot be identified. Supporting this fact, the accuracy of the vertical angle is checked before and after applying the advanced model and realized to remain unchanged ($18”$ for the Leica ScanStation P50 and $24”$ Maptek I-Site 8820, with a slight sub-arcsecond improvement). However, the sensitivity between the two other horizontal gradients of refractive index $\left[\begin{array}{cc}{\displaystyle \frac{dn}{dx}}& {\displaystyle \frac{dn}{dy}}\end{array}\right]$and range refractive index $n$ for the points located within those limited fields-of-view is partially resolved ($2\%$ and$6\%$ accuracy improvement in X-coordinates and Y-coordinates, respectively (Table 5)).

Generally, points observed at shallow vertical viewing angles exhibit reduced sensitivity to refraction effects. However, the absence of noticeable improvement in the Z-coordinates offers a stronger sensitivity between the vertical refractive index gradient ${\displaystyle \frac{dn}{dz}}$and the range refractive index $n$, which becomes more prominent at steeper vertical viewing angles (Table 6). This behavior reflects the directional dependence of 3D point coordinate accuracy on the vertical gradients of the refractive index, depending on the vertical observed angle. Conventionally, the impact of atmospheric refraction on the X- and Y-coordinates is minimal for the points near the zenith, whereas Z-coordinates are more influenced by refraction changes (i.e., $sinv$ in Equation 1).

Results underscore that shorter baselines combined with larger vertical angles (dataset 2: dam site) are more influenced by atmospheric distortions than longer baselines integrated with smaller vertical angles (dataset 1: mine site). It is implied that the improved correction model not only enhances the accuracy of the vertical angle (Table 6), but it also upgrades the range accuracy (due to diminished sensitivities (Table 7)). Those mutual effects exceed the overall 3D point coordinate accuracy. Accordingly, the advanced physical model expands corrections along 3D point coordinates consistently, where higher percentage of improvement is expected for X- and Z-direction than Y-direction. It also acknowledges that the impact of the horizontal gradients of the refractive index on the overall accuracy of 3D point coordinates is minimum, in comparison with the maximum impact of the vertical gradient of the refractive index on 3D point coordinates.

In summary, the findings indicate a moderate level of improvement with the advanced physical refractive index model, Ciddor’s model (Equations 19 and 24), where 3D point coordinate accuracies were enhanced from the centimeter to the millimeter level for the Leica ScanStation P50, and predominantly in the Y-coordinate for the Maptek I-Site 8820. However, two major concerns remain and are worth further investigation: (1) the limited parameterization of physical refractive index models for practical long-range terrestrial laser scanning (whether the Ciddor or the Closed Formula), and (2) the potential unreliability of 3D spatial gradients of atmospheric parameters, particularly vertical temperature gradients across stratified atmospheric layers, when applied to estimate refractive index gradients along the entire path in the advanced physical refractive index model. These limitations have likely contributed to the insufficient improvements reported in Table 7. To overcome these issues, hybrid physical–data-driven neural network models are proposed, supported by validated outcomes from physical modeling and cross-referenced with the control network.

5.3. Hybrid Refractive Index Model

The argued insights into physical refractive index modelling highlighted that a moderate reduction in the uncertainty of 3D point coordinates can be optimized, provided that the maximum possible precision for the estimation of the refractive index is assigned along the travelling path, with the least potential sensitivity between its spatial gradients.

Following this principle, a hybrid refractive index model is proposed. This approach enables the model to assign variable weights to the 3D spatial gradients of the refractive index, performing as a comprehensive solution to compensate for the inherent limitations of physical algorithms and to accurately evaluate atmospheric refraction along the line of sight. Consequently, the results aid in validating the field results obtained from advanced physical models and improving the overall 3D point coordinate accuracy (i.e., consistent millimeter level relative precision for the 3D point coordinates) through the nonlinear treatment of refractivity along the path across different scanning environments and using various scanners.

To support accurate prediction through this data-driven approach, it is recommended that remaining systematic errors be eliminated from both field datasets. However, they are expected to be insignificant in comparison to the priori standard deviation of observations – referring to Table 1. Generally, the following steps must be taken:

• Generate vectors of in-situ atmospheric recordings (e.g., air temperature, atmospheric pressure, and/or relative humidity and their spatial gradients) and intrinsic scanner characteristics (wavelength number, range, and angular accuracy) as the input data, to predict the refractive index as the output.
• Implement the training of a neural network for the output function of refractive index and its spatial gradients.
• Perform symbolic regression on the neural network outputs using physical interpretation of basis functions, such as Ciddor (Equations 13–19, 24, and 30) [22,23].
• Derive closed-form symbolic expressions for refractive index, applicable to new atmospheric input conditions (optional) (For illustration, the MATLAB implementation of these steps is presented in the Appendix D.).

Through symbolic expressions derived from the neural network, the increased precision for the estimation of refractive index and its gradients further reduces the sensitivity to spatial gradients of the refractive index. The posteriori uncertainties in range from the mine site dataset and vertical angle from the dam site dataset were reduced to $1.8mm$ and $9″$ for the Leica ScanStation P50, and $2.3mm$ and $15″$ for the Maptek I-Site 8820. Additionally, those consistently diminish the uncertainty of 3D point coordinates to reliable millimeter level (Figure 9).

Figure 9. Accuracy comparison of 3D point coordinates under the three proposed methods.

Figure 9. Accuracy comparison of 3D point coordinates under the three proposed methods.

Figure 9 illustrates the improvement in 3D point coordinate accuracy obtained through hybrid atmospheric correction models based on Ciddor developments - advancing refractive index estimation precision along the laser path and reducing uncertainty from the centimeter to the reliable millimeter level for the 3D point coordinates consistently. Worth emphasizing, both models, the advanced physical model and the hybrid model, consistently enhance 3D point coordinates accuracy, with particularly notable improvements in the X- and Z-coordinates. These results confirm the effectiveness of the proposed models in mitigating atmospheric errors in long-range terrestrial laser scanning compared to conventional physical refractive index modelling.

5.4. Discussions on Atmospheric Refraction Corrections

Previously, the impact of atmospheric refraction on the 3D point coordinates was investigated at two monitoring sites in terms of range and vertical angle refraction variations. Given Table 5 and Table 7, it was shown that the improved atmospheric correction model primarily affects the Z-coordinates, with the maximum impact for the points close to the zenith, while the X- and Y-coordinates are minimally influenced for points close to the horizon and negligible changes close to the zenith. In the discussion section, the comparison of all observed points in each dataset with their corresponding ranges, vertical angles, and Z-residuals - obtained through three proposed refractive index models - is presented (Figure 10 and Figure 11). This further enables visualization of the behavior of the entire 3D point cloud after applying the atmospheric correction models for the dam site dataset, acquired by the Leica ScanStation P50 with extended simulated scanning ranges.

The results from the mine site demonstrate that the hybrid physical model substantially outperforms both the conventional and advanced physical models. As outlined before, the required corrections are at a sub-millimeter level for the observed ranges. Furthermore, the conventional model shows significant variations in atmospheric corrections, particularly at longer baselines than $200m$, while the advanced model reduces overall systematic bias for the entire range datasets consistently. Note, the larger variation in conventional methods is caused as the result of linear assumption of refractive index modelling which was addressed by weighted gradient indices in advanced and hybrid methods.

The analysis of the second dataset also reveals that the observations at larger vertical angles (particularly greater than $\approx 60°)$ have a considerably higher impact than the observations at shallower vertical angles on the overall accuracy of the 3D point coordinates. This outcome arises due to the increasing sensitivity between range refraction and vertical gradient refraction, representing that the vertical gradient of refractive index dominates long-range TLS atmospheric error modelling, compared to the range refractive index in the mine site dataset. Importantly, the requirement to account for non-uniform atmospheric stratification and the anisotropic distribution of noise at larger vertical viewing angles is essential. Therefore, based on the successful results accomplished from advanced refractive index model, a combination of optimized weightings using the vertical gradient of refractive index into the advanced physical model provides a robust and intuitive solution to mitigate systematic atmospheric errors, allowing reliable millimeter- to sub-millimeter-level precision in 3D point cloud reconstruction for TLS-based deformation monitoring applications (Figure 11).

Subsequently, Figure 12 shows that, since atmospheric corrections are more significant at larger vertical angles, the high proportion of range corrections occurs from the Z-direction rather than the X- or Y-directions (i.e., the vertical angle observations are more sensitive to refractive index variations along the laser path). Consequently, this makes vertical angle corrections critical for ensuring the geometric accuracy of 3D point coordinates at long ranges (Figure 13). To further illustrate the practical significance of atmospheric effects, comprehensive 3D point cloud simulations are conducted for the dam site dataset across three representative scanning ranges ($200m$, $400m$, and $1000m$ (the maximum reported scanning range for the Leica ScanStation P50)), covering a vertical angle variation from $-10°$ to $81°$. The results of these simulations support the residual patterns identified in Figure 10 and Figure 11 (Figure 13). At higher vertical angles (between $60°$ and $80°$) and ranges greater than $200m$, the atmospheric correction increases by more than a factor of two – from $6.6mm$ to $15.2mm$. At the simulated maximum range of $1000m$, this correction reaches $51mm$.

6. Conclusions

The current research aimed to investigate the effect of atmospheric variations along the line of sight when long-range terrestrial laser scanning is required. The traveling laser beam through the atmosphere is predominantly affected by several optical phenomena. One of the most significant occurrences is the refraction of the optical path. Refraction causes deviations from the theoretical optical path and introduces variations in the intersection between the laser beam and the surface of the targets. This is identified as one of the systematic error sources in long-range scanning. To address the refractivity patterns along the optical path, the detailed mathematical developments of two physical refractive index models (the Ciddor and the Closed Formula) are elaborated here.

The physical model establishes the relationship between atmospheric conditions (air temperature, atmospheric pressure, and relative humidity) and the refractive index using the given wavelength. In the conventional approach, the average of refractive indices, based on the first and last terminals of the sightline, is obtained to address the linear condition of refractivity along the path. Furthermore, a higher precision on the order of $\left[1-5\right]\times {10}^{-8}$is recommended for the estimation of refractive index to ensure millimeter or sub-millimeter accuracy for the observations [1]. For advancements of the physical models, 3D spatial gradients of the refractive index are developed in the current work, in relation to the nonlinear physical modelling where varying refractive indices along the path are demanded. These requirements are applied when the laser pulses experience the different medium (vertically stratified atmospheric layers).

In the field experiments for the atmospheric modelling, two long-range terrestrial laser scanners (Leica ScanStation P50 and Maptek I-Site 8820) were employed for quality testing under varying atmospheric conditions at two deformation monitoring sites: a mine site - enabling the long-range scanning over $800m$ - and the dam site - providing the flexibility of steep vertical viewing angle close to zenith $80°.$ Both sites were controlled under a calibrated network arrangement, utilizing an on-site range consistency calibration method supported by GPS control points, along with field survey observations to verify vertical angle accuracy. During each scanning setup, the atmospheric conditions of air with respect to individual scanner station (in-situ atmospheric recordings attached to each station) were tracked across the entire sites.

The recognized ranges from the mine site dataset were evaluated for network reliability, and after maximizing the precision of refractive index estimation through advanced approach of physical refractive index model, the improvements of$34\%$ for the Leica ScanStation P50 and 16% the Maptek I-Site 8820 in range accuracy were observed. Since a number of points were located at the shallow vertical angle in the mine site, the sensitivity to 3D spatial gradients of the refractive index – specifically, vertical gradient of refractive index across stratified atmospheric layers – was not detected, and this led to limited improvement in 3D point coordinate uncertainties. The dam site dataset demonstrates mitigation strategies for this sensitivity, particularly for points located at higher elevations (with vertical viewing angles larger than 60°). With the advanced physical model, improvements occurred not only in the a posteriori accuracy of the vertical gradient of the refractive index (44% (from $18”$ to $10”$) for the Leica ScanStation P50 and $20\%$ (from $24”$ to 19”) for the Maptek I-Site 8820), but also in the accuracies of the 3D point coordinates (reliable centimeter- to millimeter-level accuracy). The reason is that the points, located on approximately the same horizontal plane as the scanner station, are more uniformly affected by atmospheric corrections (dataset 1: mine site), while the points, located at different horizontal planes, require larger X- and Z-coordinate corrections (dataset 2: dam site). The hybrid model, referencing back to the advanced physical refractive index model, is advised to achieve higher millimeter-level relative precision in 3D point coordinates consistently. Millimeter-level accuracy is ultimately attained by reducing the range and vertical angle uncertainties to $1.8mm$ and $9″$ for the Leica ScanStation P50, and $2.3mm$ and $15″$ for the Maptek I-Site 8820 (Table 6). The advantage of the proposed approaches - the advanced and hybrid physical models - is their ability to represent real-world refractivity conditions along the laser path by maximizing the precision of refractive index estimation and minimizing sensitivity to spatial gradients, achieved through a stochastic weighting of the vertical refractive index gradient within different atmospheric layers.

For future work, several pathways exist toward achieving millimeter- or sub-millimeter-level accuracy in the field calibration of long-range terrestrial laser scanning. First, additional optical effects of the laser line, such as scattering and reflection, play an integral role in determining robust radiometric and spatial calibration results. Second, the algorithmic steps presented here for the physical refractive index model are recommended to be further modified based on more accurate in-situ atmospheric observations (e.g., temperature accuracy better than $0.5°C$). This is particularly suggested for highly sensitive and high-risk deformation projects. Preferably, attached thermometer arrangements should enable recording of the epoch-wise atmospheric conditions of the laser. Alternatively, for advanced physical refractive index modelling, it is strongly recommended to create a temperature heat map sensitive to height profiles, as the dominating factor, to better indicate temperature variations across the monitoring test sites. Third, the use of an appropriate stochastic model for 3D spherical observations is highly advised for comprehensive system calibration of the 3D point cloud, addressing geometric error models. Finally, the relevant atmospheric correction factors should be ideally applied within the instrument, enhanced by a robust atmospheric measurement technique along the entire path which can be supported by the manufacturing principal assembly.