Image shift due to atmospheric refraction: prediction by numerical weather modeling and machine learning

Abstract. We develop and study two approaches for the prediction of optical refraction effects in the lower atmosphere. Refraction can cause apparent displacement or distortion of targets when viewed by imaging systems or produce steering when propagating laser beams. Low-cost, time-lapse camera systems were deployed at two locations in New Mexico to measure image displacements of mountain ridge targets due to atmospheric refraction as a function of time. Measurements for selected days were compared with image displacement predictions provided by (1) a ray-tracing evaluation of numerical weather prediction data and (2) a machine learning algorithm with measured meteorological values as inputs. The model approaches are described and the target displacement prediction results for both were found to be consistent with the field imagery in overall amplitude and phase. However, short time variations in the experimental results were not captured by the predictions where sampling limitations and uncaptured localized events were factors.


Introduction
The Earth's atmosphere includes several phenomena that affect the propagation of light. While scattering and absorption by clouds, fogs, and aerosols primarily affect the intensity of the radiation received, the atmosphere can also affect the spatial resolution properties and the propagation trajectory of light. For example, atmospheric turbulence causes image shimmering and blurring and is stochastic in nature with fluctuations over short timescales (e.g., milliseconds). Another phenomenon is atmospheric refraction where refractive index gradients can steer or bend light rays. The index gradients are associated with changes in air density, which for optical wavelengths is primarily a function of air temperature gradients. Atmospheric refraction tends to cause more deterministic, larger-scale effects than turbulence and the effects can persist from minutes to hours. [1][2][3][4] The interest here is refraction in the lower atmosphere, which can cause apparent displacement or distortion of objects when viewed by imaging systems or produce steering when propagating laser beams.
For several years, we have been developing a low-cost, mobile camera system to study atmospheric refraction at New Mexico State University. One system was recently deployed at White Sands Missile Range (WSMR), New Mexico (NM), and a second system was set up at the Jornada Experimental Range (JER), near Las Cruces, NM. Both systems collect time-lapse images of distant natural targets, such as mountain ridges. A time-lapse system was previously used in Las Cruces, New Mexico, with a building as a target to study diurnal image displacement due to refraction. 5 A similar system in Dayton, Ohio, was used by Basu et al. 4 to investigate the temporal variations of the refractive index gradient. Time-lapse imagery has also been used to investigate the apparent stretch and compression of objects due to atmospheric refraction lensing effects 6 and the approach has also been applied to the estimation of turbulence strength. 7 The prediction of atmospheric refraction effects can be advantageous for many terrestrial optical applications where prior knowledge of the light's trajectory can improve the speed and accuracy of pointing and tracking functions. The goal of this paper is to develop and evaluate two different methods, numerical weather prediction (NWP) and machine learning (ML), for predicting image displacement due to atmospheric refraction. NWP is an attractive approach for our application as it is deeply rooted in physics. However, it is computationally expensive, and the results are subject to initial conditions and terrain characteristics. An alternative, more empirical tactic is to apply an ML algorithm to build a predictive model based on local meteorological data. In this paper, we describe our application of NWP and ML methods to image displacement due to refraction and compare the results with time-lapse camera measurements.

Time-Lapse Image Collection and Processing
During January and February of 2018, we collected image data with a time-lapse camera located at WSMR that was pointed generally north at a natural desert landscape and a mountain range (Oscura Mountains) on the horizon. Another camera was set up at JER and was pointed west to image a mountain range (Dona Ana Mountains) and a desert valley. This system began collecting images in May 2018 and is still operating. The mountain ridgelines observed were at distances of about 20 km for JER and over 100 km for WSMR.
The battery-powered camera systems are easily transportable and consist of a weatherproof case on a tripod that contains a Nikon D5200 camera operated in a time-lapse mode. A zoom lens is set at its maximum focal length of 400 mm for the WSMR measurements and at 300 mm for the JER observations. The camera is typically programmed to collect images in 5-min intervals with a fixed 5.6 f-number and automatic shutter speed. Example frames of the mountain targets and valleys for the WSMR and JER experiments are shown in Fig. 1. The rectangles indicate areas in the images that were cropped and used for the refraction analysis in this paper.
Local weather has a significant effect on the vertical temperature gradients that are primarily responsible for the atmospheric refraction effects. The weather variables of interest in our study include temperature, humidity, pressure, and solar radiation. For the WSMR experiment, online metrological data were downloaded from a weather station near the target mountain. A Davis Vantage Vue brand weather station next to the camera was utilized for the JER experiment. These measurements are interpolated in time to align with the time-lapse image frames. The Kanade-Lucas-Tomasi point-tracking algorithm 8,9 is implemented to measure the apparent motion of the mountain ridges in the images. 10 The general steps of the data-processing approach are depicted in Fig. 2. An area containing the far-field target in each frame is cropped and the N "best" features in the subframe are determined by a threshold setting. The results are stored in a feature list in descending order of "goodness." Then, the algorithm tracks these features in consecutive frames. A near-field object reference close to the camera is also selected and point-tracking of this feature is applied in the analysis to remove shifts in the far-field images that are due to camera platform motion. Figure 3(a) shows some selected points associated with the cropped image of the mountain peak in the JER experiment, and Fig. 3(b) shows the positions of these points in the next consecutive frame. The point-tracking algorithm on average selects the same points in each frame. The vertical positions of the multiple points are averaged to give a measurement of the ridge position.
After point-tracking, the near-field average pixel coordinates are subtracted from the far-field coordinates frame by frame to obtain the apparent position of the far-field target. Displacement of the target's apparent position from frame to frame is attributed to changes in atmospheric refraction. The most significant shifts are found to occur in the vertical direction.

Numerical Weather Prediction and Ray Tracing
NWP is a discipline where governing equations and parameterizations that describe fluid flow and other physical processes are applied to current (or previous) weather conditions to provide a future forecast. For our purposes, the model results can be used to predict the vertical structure of the refractive index in the atmosphere. However, an extension of the established models is required to provide higher spatial resolution along our paths of interest. 11,12 In this section, we describe our approach for using refractive index gradient data generated by NWP and the application of ray tracing to determine corresponding image shifts. The results of the approach are compared with time-lapse measurements in Sec. 3. The numerical weather model (called the Weather Research and Forecasting -WRF model) uses initial and boundary conditions from a global-scale reanalysis data set (called ERA-5) and topographical effects to generate the refractive gradient data for a particular location and time range corresponding to our field measurements. 12 Figure 4 illustrates the refractive index gradient ½dnðhÞ∕dh model results for the WSMR experiment time-lapse imaging path on the morning of February 5, 2018, where nðhÞ is the vertical profile of refraction index as a function of height h. The gradient values are presented as a function of altitude relative to mean sea level and a function of distance along the propagation path. In Fig. 4, the camera site is on the left and the mountain ridge is on the right.   The camera site is on the left and the mountain ridge is to the right. The target corresponding to this result is the rectangular area indicated at the lower right side of the mountain ridge in Fig. 1(b). In fact, the actual mountain peak in Fig. 1(b) is a narrow ridge that is not adequately resolved by the current NWP model spatial resolution (∼1 km). Thus, NWP results are not available for the mountain peak. However, the model result shown in Fig. 5 allows us to examine the shift of a lower portion of the mountain ridge at about 1400 m in elevation where the viewing path is nearly horizontal across the basin.
Ray tracing through the gradient profiles is used to determine the image displacements predicted by the model. Ray-tracing techniques are often applied for refraction analysis over near-ground horizontal paths assuming diffraction effects are not significant. 3,13 The NWP data essentially consist of "blocks" of constant refractive index gradient, as indicated in Figs. 4 and 5. Rather than using a conventional linear ray-tracing algorithm that requires subsampling the blocks to provide accurate trajectories, we apply a second-order ray tracer 14 where the linear ray transfer equation is expanded with a quadratic correction term to model the curved ray trajectory within each block. This approach requires only one tracing step for each data block and is significantly faster than a linear ray trace approach with subsampling. A summary of this method is now presented.
Consider a two-dimensional form of the eikonal equation describing the ray trajectory in an inhomogeneous media: 15   (1) where h is the height and x is the horizontal distance. This expression assumes horizontal paraxial propagation, and the refractive index gradient is vertical. Each block in the NWP data has a constant vertical gradient value indicated by κ ¼ dn∕dh and, assuming nðhÞ ≈ 1, the solution to Eq. (1) for the ray vertical position within a block can be given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 5 1 hðxÞ ¼ where θ 0 and h 0 are the initial ray angle and height, respectively. Taking the derivative of Eq. (2) with respect to x gives the ray trajectory angle as a function of distance, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 8 4 With respect to the block boundaries in the NWP index gradient data, Eqs. (2) and (3) are used in succession to transfer the ray height between one boundary and the next and then provide the bending of the ray angle for the next block. Iteratively, the equations become: and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 4 6 1 where Δx is the distance between adjacent blocks (∼1 km for our NWP data) and j is an index that identifies the different blocks. Example ray trajectories generated by the tracing approach are illustrated in Figs. 4 and 5. Rays (200 rays for WSMR and 100 rays for JER) are launched from the image target (mountain ridge elevation of ∼2200 m for WSMR and the basin edge elevation of ∼1415 m for JER) for a range of initial angles (−0.1 to −8 mrad for WSMR and −1 to þ1 mrad for JER). The ray trajectories are traced through the model gradients until they reach the ground near the camera. The specific ray that strikes the ground at the camera location is identified and a line is backprojected at the ray arrival angle. The height of this line at the target plane indicates the apparent target position as seen by the camera. The apparent positions are calculated for successive model frames and the relative shifts are determined. For the results presented here, NWP results were computed for 10-min intervals, and the ray-tracing procedure was applied. We note that it is also possible to trace the rays from the camera position toward the target position.

Machine Learning Predictions
In this section, we describe an ML approach to predict image displacement due to atmospheric refraction based on a set of measured metrological values. The input variables we use for prediction are temperature (T), humidity (H), pressure (P), and solar radiation (S). The predicted output is the image displacement (ŷ) due to refraction. Other available meteorological data, such as wind speed, were also applied as test inputs to the ML model, but we found that these alternative parameters had little influence on the prediction results. The weather station for the JER experiment provides measurements of these variables at 15-min intervals. In addition, we utilize other local measurements available online at hourly intervals. Prior to input to the algorithm, the measurement values are normalized to the range of (0,1) by dividing each value by the maximum value observed. Our ML prediction process follows the conventional approach of splitting the image displacement and meteorological data into three sets: training, validation, and testing. Our prediction is a comparison of the ML results with the testing data set.
The ML approach is based on linear regression and we assumed a model of the form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 7 2 3ŷ ðT; H; P; S; wÞ ¼ w 1 þ w 2 T þ w 3 T 2 þ w 4 H þ w 5 H 2 þ w 6 P þ w 7 P 2 þ w 8 S þ w 9 S 2 þ w 10 TH þ w 11 TP þ w 12 TS þ w 13 HP þ w 14 HS þ w 15 PS; (6) where w ¼ ½w 1 : : : : : : : : : : : : : : : w 15 T are the coefficient weights, and linear, square, and pairwise products of the meteorological parameters are used as nonlinear kernel functions. The choice of kernel functions was based on trial-and-error experimentation. The coefficient values w are determined by fitting Eq. (6) to the training data by minimizing an error function that measures the misfit betweenŷ and the measured values y ¼ ½y 1 : : : : : : : : : : : : y N T as a function of w. Our choice of error function is the regularized squared error: 16,17 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 0 0 ½ŷðT n ; H n ; P n ; S n ; wÞ − y n 2 þ λ 2 kwk 2 : The term ðλ∕2Þkwk 2 is a penalty (regularization) term to avoid overfitting where the parameter λ is an input to the model that governs the relative importance of the regularization term compared with the squared error term. Given N data points ðT n ; H n ; P n ; S n ; y n Þ, the coefficients w that minimize the cost function in Eq. (7) are obtained in closed form by differentiating EðwÞ with respect to w, setting the result to zero, and solving for w. This produces the following well-known result: 17 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 4 7 8 where I is the identity matrix and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 4 3 4 As shown in Eq. (9), the input observations are arranged as row vectors. The algorithm steps involve constructing the matrix Φ, obtaining the vector y; and applying Eq. (8) to compute the weights w. Referring to Eq. (7), in order to generate a predictive model that generalizes well for new data input, the selection of the value of λ needs to strike a balance between overfitting and underfitting of the training data. A larger value of λ tends toward a "simpler" fit but there is more likelihood of underfitting the data where some of the dominant trends are not captured. On the other hand, a low value of λ provides a more "complex" fit but there is more likelihood of overfitting the data where specific noise events are captured as trends. We use a tuning approach to determine the value of λ where a search over a range of values is performed. For each λ trial, the model is first fit to the training data set. This model realization is next applied to the validation data, and finally the mean squared error (MSE) between the model and measured validation data displacement result is calculated. The value of λ that provides the lowest MSE is selected and used for the next step of prediction comparisons using the testing data.
Finally, we mention that for the JER experiment, the model form in Eq. (6) was used but we added a fifth binary variable (D) that takes either a value of 1 if the sky is "clear" or 0 if it is "cloudy." The value of D is determined from a visual analysis of the sky conditions in the timelapse imagery. If the sky appears to be cloudy in more than 50% of the frames during the daytime hours, then D is set to 0, otherwise it is set to 1. Although this is an unrefined measurement and a simple binary parameter, we found the approach improved the accuracy of the model for varying sky conditions.

Results
NWP results were computed for 10-min intervals and the ray-tracing procedure was applied to determine the target shifts. The shift results were interpolated in time to align with the time-lapse image frames. Figure 6 shows comparative results for the WSMR experiment on February 5, 2018, as a function of time-of-day where the red curve is the apparent shift (in radians) of the mountain ridge as derived from ray tracing the NWP data. The black curve is the predicted shift by the ML algorithm where, prior to the date shown, we trained the algorithm on 4 days of displacement and meteorological data (500 data points) and tuned the result on 2 days of validation data. The search range for λ was (0, 40) and, for these results, λ was found to be 2. The blue curve is the shift measured in the actual camera frames. Measurements are only available for the daytime hours as the ridge is not clearly distinguishable by the camera at night.
The general downward drift throughout the daytime hours, as shown in Fig. 6, is an effect we commonly observe in clear weather and this corresponds to a slow reduction in the average refractive index gradient of the atmosphere along our line of sight. The NWP ray trace and the ML model results agree well in amplitude and phase with the image shift results from the camera measurements, although there are differences in the short-time variations. By manually evaluating some of the time-lapse images, we verified that the point tracking results appear to be accurate, and so we believe the short-time variations are due to the atmosphere. These excursions probably represent turbulent fluctuations in the refractive index that are not completely captured by the NWP and it is not possible to predict these short time variations by the ML model due to the differences in sampling where the time intervals for the meteorological data (1 h) are much longer than the time-lapse intervals (5 min). The meteorological measurements were also not collected directly in the imaging path, which could contribute to further differences in the measurements and ML results. Figure 7 presents example shift results for the lower part of the ridge in the JER experiment on July 18, 2018. Details of the data-processing approach are the same as described for the WSMR results (Fig. 6). For this result, the training and validation data sets consisted of 2 days each and λ was found to be 0.01. Like the WSMR results, the NWP ray trace and ML predictions in this case demonstrate the same overall amplitude and phase behavior as the time-lapse measurements, but the short time variations are not captured. Figure 8 shows the measurement and ML prediction results for the mountain peak in the JER experiment [ Fig. 1(b)] over a period of 6 days from January 20, 2019 to January 25, 2019. As discussed in Sec. 2.2, NWP results are not available for this peak target. The 4-day training set in this case included different weather conditions (sunny, very cloudy, and mixtures of sunny and cloudy conditions). The value of λ was found to be 5. The weather was sunny and clear for all the days shown except for January 22 when the sky was cloudy. The results show that the ML model prediction for different weather conditions and over different path angles is consistent with the general trends of the measurements. We note that the cloudy-day result shows less of the daily downward progression of the peak position. Note that the overall shifts for these JER mountain peak results are significantly smaller than the results for WSMR (Fig. 6). This is likely because the WSMR ridgeline is much further from the camera (>100 km) than the ridgeline for JER (∼20 km).
We end this section with a few comments about the range of values found for the regularization parameter λ. Generally, we expect a smaller value of λ (less regularization) to be found  when the short time variation amplitudes are relatively smaller than the mean trend deviations within the training and validation data. However, the value of λ selected through our tuning approach appears to be sensitive to other factors such as the shape of the mean trend curve. Also, because we trained and validated the model on relatively small data sets, even adding a few days of data or including a few unusual shift results can affect the λ value. It is important to note that although the tuning process identifies a particular λ value, we found through additional testing that there is a range that typically produces nearly the same prediction result. For example, similar ML results to that shown in Fig. 8 can be produced with λ values ranging from 0.001 to about 10.

Conclusions
We found that NWP along with a ray-tracing technique can be used for prediction of image displacement effects due to atmospheric refraction. The ray-tracing approach that we applied to determine the effect of the gradients produced by the model was straightforward to implement and provided credible results. The model results of target displacement were consistent with field imagery in amplitude and overall trend. However, the model could not predict some of the short time variations in the field measurements, which may be due to localized events that are not completely captured by the numerical model grid or simulation process. As an alternative to NWP, we explored the use of an ML algorithm to build a predictive model based on meteorological data collected near the camera location. We found that our ML model was successful in making predictions over different weather conditions. Similar to the NWP result, the ML model prediction could not follow the short time excursions of the field results. In this case, the slow sample rate of the meteorological data compared to the time-lapse image frame rate is an aggravating factor. We are now working to determine if the ML approach can be extended to encompass different seasons as well as different weather conditions. We also are investigating our ability to apply image point-tracking and ML prediction to detect geometrical distortions of the target image rather than just the apparent target shift.