Real-time sea-level monitoring using Kalman filtering of GNSS-R data

Current GNSS-R (GNSS reflectometry) techniques for sea surface measurements require data collection over longer periods, limiting their usability for real-time applications. In this work, we present a new, alternative GNSS-R approach based on the unscented Kalman filter and the so-called inverse modeling approach. The new method makes use of a mathematical description that relates SNR (signal-to-noise ratio) variations to multipath effects and uses a B-spline formalism to obtain time series of reflector height. The presented algorithm can provide results in real time with a precision that is significantly better than spectral inversion methods and almost comparable to results from inverse modeling in post-processing mode. To verify the performance, the method has been tested at station GTGU at the Onsala Space Observatory, Sweden, and at the station SPBY in Spring Bay, Australia. The RMS (root mean square) error with respect to nearby tide gauge data was found to be 2.0 cm at GTGU and 4.8 cm at SPBY when evaluating the output corresponding to real-time analysis. The method can also be applied in post-processing, resulting in RMS errors of 1.5 cm and 3.3 cm for GTGU and SPBY, respectively. Finally, based on SNR data from GTGU, it is also shown that the Kalman filter approach is able to detect the presence of sea ice with a higher temporal resolution than the previous methods and traditional remote sensing techniques which monitor ice in coastal regions.


Introduction
A significant portion of the world's population lives and operates in coastal regions and is thus susceptible to hazards originating from the state of the sea (Neumann et al. 2015). To detect imminent threats such as storm surges and tsunamis, real-time measurements of the current sea level are necessary (Holgate et al. 2008). Real-time sea-level reports are also necessary for day-to-day operations of, for example, shipping routes and ports (Pugh 2004). For these reasons, the Global Sea Level Observing System (GLOSS) initiative of the Intergovernmental Oceanographic Commission (IOC) of UNESCO encourages the collection and distribution of data in real time or near real time (UNESCO/IOC 2012). In the implementation plan of GLOSS, the IOC also calls for the usage of GNSS stations to monitor the movement of tide gauges. Since it has previously been proved that the GNSS stations can be used directly to monitor sea level ) with a technique called GNSS reflectometry (GNSS-R), we look into using this technique to retrieve sea level in real time.
The fundamental idea of GNSS reflectometry relies on the principle of using multipath effects to deduce properties about the surroundings of a GNSS antenna. Multipath effects are caused by a variety of properties of the environment around a GNSS antenna (Nievinski and Larson 2014a). As such, the multipath effects, which are considered a major error source for positioning and navigation applications, also contain some information about their causes. GNSS reflectometry works like a passive radar system in that it uses signals that are broadcasted by some other source, in this case the GNSS satellites, and utilize these for measurement purposes other than those intended originally. All GNSS signals can be used for this purpose, and with more GNSS systems and satellites still being deployed, GNSS reflectometry is, therefore, becoming increasingly more interesting for remote sensing and related fields. Reflected GNSS signals can be recorded both from airborne or space-borne platforms (Martin-Neira 1993) and from ground-based systems (Anderson 2000). However, their principle, as well as their purpose, is quite different, and the focus of this work will be the latter, i.e., installations near the ground. More specifically, the focus is set to the interferometric approach to GNSS reflectometry, since this technique uses only one single off-the-shelf receiver and antenna pair (Bilich and Larson 2007), keeping costs relatively low and giving access to data from already existing networks of GNSS installations (Dow et al. 2009).
The most common method for using GNSS signals for GNSS-R is through analyzing their signal-to-noise ratio (SNR) with spectral methods (Jin et al. 2014), e.g., with Lomb-Scargle analysis, which can reveal properties such as the snow depth around the antenna (Larson et al. 2009) or the soil moisture of the surrounding ground (Larson et al. 2008). However, the method has been proven less suitable when, for example, the reflection surface is dynamic, as is the case for sea-level measurements, at which point a correction term is needed .
As an alternative approach for measuring on dynamic surfaces, we have previously developed an inversion-based method (Strandberg et al. 2016) which is also using SNR from off-the-shelf equipment and existing GNSS installations on the ground. This method builds on the idea of fitting an inverse model that represents the recorded SNR data. Thus, by using least-squares methods, it is possible to retrieve properties of the reflecting surfaces from the SNR data. We have shown that this method retrieves the sea surface height with better precision than Lomb-Scargle methods (Strandberg et al. 2016). However, a drawback of the inverse modeling method with least-squares analysis is that it requires that all data for the period of interest are available when the inversion starts. This can create significant delays between the time when data are recorded and when the final analysis results are presented. In cases, when data in near real time are required, the method will be inadequate.
Because of this drawback, we introduce a new method to retrieve information from GNSS SNR data, which addresses the aforementioned latency issue by using Kalman filtering techniques together with the inverse modeling formalism. With this new method, it is possible to retrieve precise results in real time, retaining the detailed functional description of the inverse modeling, and therefore opening up for new applications in the field of GNSS reflectometry.

Current state of GNSS reflectometry
GNSS reflectometry is based on the idea of accessing the sources of multipath through their impact on the SNR of the GNSS signal. As the coherent part of the reflected signal interferes with the direct signal, the resulting SNR varies in strength according to a characteristic pattern that depends on the elevation of the GNSS satellite and other parameters (Nievinski and Larson 2014a).
Under the assumption of a locally flat and horizontal reflector, the SNR as a function of elevation can be expressed as a monotonically increasing trend superimposed by a high-frequency oscillation pattern (Bilich and Larson 2007), as presented in Fig. 1. After removing the trend with a low-order polynomial, the high-frequency residual oscillations behave like a damped sinusoid (Strandberg et al. 2016), which can be described by where is the elevation to the GNSS satellite corrected for atmospheric refraction (Bennett 1982), k the wave number of the signal, and h the vertical distance between the antenna and the horizontal reflector. The remaining parameters, A , , and , are aggregate variables depending on several properties, including, for example, antenna gain patterns of the receiver installation and the satellite, and electromagnetic properties of the reflector (Nievinski and Larson 2014a).
Under the assumption of a static reflector, we can see that the term 2kh relates to an oscillation frequency of the sinusoidal part of the SNR, depending only on sin . Therefore, since we know the elevation angles of the satellites with sufficient precision from satellite broadcast orbits, we can use the oscillation frequency of the SNR from a single pass of a GNSS satellite to determine the (1) SNR = A sin (2kh sin + ) ⋅ exp 4 k 2 sin 2 , reflector height h . Depending on the height of the antenna above the reflecting surfaces, roughly 20 min of SNR data are enough for a single height retrieval (Williams and Nievinski 2017). This method, which we will refer to as spectral analysis, has been successfully demonstrated to be able to retrieve, for example, snow depth , where the assumption of a static reflector is appropriate. However, when the reflector height is dynamic, some adjustments are needed, for example, by adding the first-order correction term for the changing sea level . Still, unless the retrieved sealevel values are used to form temporal averages, the precision does not live up to the performance requirements of on an operational tide gauge (Larson et al. 2017).
Another way of retrieving parameters from the SNR measurements is by using inverse modeling (Nievinski and Larson 2014b), for example, by fitting (1) or a similar approximation function directly to the SNR data either by interval analysis (Reinking 2016) or by least-squares adjustment (Strandberg et al. 2016). In doing so, it is possible to go from a static reflector height to a dynamically changing one by modeling temporal variations of the reflector height with, for example, a set of B-spline functions (Strandberg et al. 2016). Another improvement with the inverse modeling compared to spectral methods is that a more realistic model of the measurements is used. Whereas the spectral analysis only accounts for the oscillation frequency of the SNR, the inverse modeling also considers the shape of the oscillations. Therefore, the inverse modeling allows retrieving more of the intrinsic information that the data carry, than spectral methods do. In contrast to using only data from a single satellite pass to determine one measurement of the reflector height as the spectral method does, the inverse modeling method is also capable of combining all concurrent data in a coherent inversion step, regardless of the satellite or the satellite system. This ability to combine data from several satellites has been shown to increase the robustness and precision of the retrievals.
While inverse modeling improves several aspects of the parameter retrieval, it also comes with a drawback. As the inversion is performed using least-squares adjustment, the method is inherently off-line. Since all data for the inversion have to be available when the inversion is started, the method has significant delays before any results can be presented. The spectral methods are closer to real time, but still require a significant portion of a satellite pass before any reflector height measurements can be retrieved. This gives both a delay but also irregular intervals between reflector height estimates. Therefore, we focus on alternative methods that allow us to deliver results in real time, and in this pursuit, we make use of the concept of Kalman filtering.

Dynamic inverse modeling of SNR data with a Kalman filter
Dynamic parameter estimation in real-time applications usually relies on sequential filtering techniques. These allow updating the state variables describing the underlying physical properties using measurements that relate to the state variables either directly or through known functional expressions. Since the introduction of Kalman filtering techniques (Kalman 1960), they have been prevalent in the field of dynamic parameter estimation. The original definition of the Kalman filter built upon the assumption of linear systems for predicting and updating state estimates, making them less ideal for systems with nonlinearities (Jazwinski 1970). Therefore, as our implementation builds on the highly nonlinear relation in (1) for its update procedure, we have to resort to alternative techniques. Fortunately, the nonlinearity in the update step can also be dealt with the unscented Kalman filter (UKF) approach (Merwe and Der 2000) in which the predicted measurements and their covariances are obtained approximately after selecting so-called sigma points around the predicted state, parsing all these points through the nonlinear function h(x).
Our implementation utilizes an ordinary, linear prediction procedure, predicting the change to both the state vector and the covariance matrix where the state vector x consist of our estimation of the reflector height and the other parameters of (1). As mentioned above, the update step will follow the procedure of the UKF approach, i.e.
1. Compute 2L + 1 sigma points based on the predicted state and its covariance L corresponds to the size of the state vector and can be computed as = 2 (L + ) − L . Following the recommendations by Merwe and Der (2000), we have chosen to set = 10 −3 and = 0 in our implementation.
(2) Utilizing this update procedure ensures that the nonlinearity of (1) is taken care of without explicitly linearizing the problem.

Parametrization, process noise, and observation noise
As denoted in (1), detrended SNR can be expressed as a function of elevation angle , reflector height h(t) , damping , amplitude A , and phase . The latter two parameters are satellite system specific and must thus be estimated separately for each GNSS, as well as for each frequency. On the other hand, reflector height and damping are common to all SNR observations, which means that the combination of several GNSS will improve the estimation of all parameters and reduce the risk of having outliers bias the target parameters.
Damping, amplitude, and phase mostly depend on slowly changing parameters in the environment around a station. Therefore, they can be modeled with a random walk process that has very small process noise. However, for the sea surface height, it is hard to find a dynamic model that works for different sites, as the instantaneous sea level can be both driven by astronomical and meteorological tides to a varying extent, depending on the site and its location. Therefore, we turn instead to the B-spline method used by Strandberg et al. (2016). This approach uses quadratic B-spline coefficients in order to represent the temporal variations of the reflector height. In general, the Kalman filter state vector can only have a finite number of parameters. Thus, if we want to be able to use a Kalman filter without being limited to a fixed time interval, we have to change the definition range of our B-splines while the filter is running.

Dynamic introduction of B-spline coefficients
If temporal changes of the reflector height h(t) are described by quadratic B-spline functions B 2 i , we only need to estimate a set of scaling coefficients i . Then, we are able to evaluate the function at any point in time, i.e.
Considering now that B-spline functions are defined on the zero-degree base function and the recursive relation it is obvious that the evaluation of (18) for any higher degree k at a certain point in time t does only require the knowledge of four node values that surround the epoch of interest. However, for evaluating h(t) at any arbitrary epoch, the knowledge of all scaling coefficients i is necessary. In case of dynamic filtering, this restriction poses a problem, as the dimension of the state vector depends on the total length of the observational time series and is not predictable unless data are processed in batches. However, as this contradicts the idea of continuous Kalman filtering, we have to use another method. One method for using B-splines in a Kalman setting, presented by Jauch et al. (2017), is based on utilizing a fixed set of parameters, and periodically changing the definition range. This concept is illustrated in Fig. 2 and makes use of the fact that B-spline base functions are compactly supported, i.e., they are only valid over a certain interval.
As stated earlier, four quadratic B-spline functions are necessary to model h(t) at any given epoch. If we consider that we only need to have those B-spline coefficients in the state vector that relate to the current update epoch, the state vector is of constant size. However, as one approaches an epoch where a new B-spline base function is necessary to cover the definition range, some actions are needed to deal with this special case. As shown in Fig. 2, the four B-spline coefficients and the corresponding covariances are continuously updated as the filter moves along in time until the point t j where the new node at t j+1 needs to be introduced and the node which is no longer required ( i−3 ) has to be dropped from the state vector. As nodes i−2 , i−1 , and i are still needed, they change position in the state vector to become the first three entries for the B-spline representation, which opens the space for the new node, to be put on the fourth location in the state vector. As depicted in the second step in Fig. 2, the initialization of the new node happens with the estimated value of the third node and a slightly larger variance, which allows that following updates will correct the Fig. 2 Dynamic updates of the state vector are necessary in order to be able to take advantage of the B-spline functional expression. The three steps illustrate how the state vector and the corresponding covariance matrix are handled, before, during, and after, an epoch where a new B-spline interval is entered. The reader is referred to the related section for an indepth description of individual stages 1 3 coefficient according to the next SNR measurements. Correlations between the new coefficients and the other three coefficients are reconstructed from the covariance matrix before the shifting takes place. Having introduced the new coefficients and considered their new validity range, the filter can continue with the normal prediction/update sequence and estimate reflector height, damping, amplitudes, and phases continuously and in close real time. As the latter parameter types are modeled as random walk processes, there is no need for special considerations, except that only those satellite-specific parameters are updated when the particular satellites are in view at a certain point in time. B-spline coefficients which are no longer in the state vector are stored on disk so that one can later obtain a timer series of h(t) or evaluate the reflector height at any given point in time.

The choice of the stochastic model
While B-spline coefficients are treated as constant parameters, i.e., having zero process noise and periodically being exchanged for new nodes with some uncertainty, the other parameters are modeled as random walk processes in order to consider that environmental changes can lead to variations of these parameters. For example, the amplitudes and the damping coefficient are related to the reflectivity and thus allow concluding on parameters like ice coverage (Strandberg et al. 2017) or vegetation (Small et al. 2010). Considering that the phase of the SNR interference pattern relates, for example, to soil moisture, it is obvious that also this parameter should be allowed to vary over time.
While the filter determines estimates of the systemdependent parameters amplitude and phase for each GNSS and frequency, they are assumed to behave statistically similarly across all systems. Thus, the process noises of all amplitudes are the same, as well as the process noises for all phase parameters. This choice was made to limit the search space in order to find functional settings.
The values have been chosen to guarantee a stable solution on 20-day long data sets used for testing purposes, without any diverging parameters. The effects of changing the respective process noise variances are presented in Fig. 3. Based on these results, the chosen model parameters used for dynamic inversion of the SNR observations are presented in Table 1. As the true observation noise for the SNR data is unknown, the observation noise is estimated per satellite system and frequency from the average residual of the last hour of SNR observations, resulting in noise variances around 100-200 V per volt squared (V/V) 2 .
To be able to run the filter in real time, it is also necessary to perform the aforementioned detrending of the SNR series in real time as well. The trend depends largely on the equipment used, and the surroundings of the antenna (Nievinski and Larson 2014a), and is therefore only slowly varying in time.
Here, we introduce two alternatives which take advantage of the stability of the trends: Either the polynomial coefficients are introduced as parameters in the state vector of the Kalman filter, or the polynomial for each recurring arc is determined through a running average of the last few passages. To not increase the number of parameters in the filter overly much, we chose to use the second alternative.

Test sites
In order to assess the precision of the Kalman filter based retrievals, we used SNR data from two test sites with different characteristics. These sites have previously been used to demonstrate the least-squares inversion algorithm (Strandberg et al. 2016). The first site is the GTGU test installation at the Onsala Space Observatory, Sweden. GTGU is a purposefully built GNSS-R installation where a standard off-the-shelf receiver and antenna pair has been mounted on a pole sticking out over the water to provide optimal viewing angles of the sea surface. The antenna oversees open water in an azimuth range that is larger than 180 degrees, see the top panel of Fig. 4. The equipment consists of a LEICA AR25 choke ring antenna paired with a LEICA GRX1200 receiver recording GPS and GLO-NASS data at a rate of 1 Hz. Close to GTGU, a stilling well tide gauge is situated a few hundred meters from the antenna. This tide gauge has been used in the validation process as a reference and is assumed to provide ground-truth values. The tide gauge logs water level with millimeter accuracy once every minute using three different types of instruments, one laser, one radar, and two bubble pressure sensors (Pugh 1972), and is described by Wahlbom (2015). The latter are susceptible to biases caused by changes in the water salt content, and the radar has experienced problems with false detections. Therefore, we decided on using the laser sensor as the reference.
Tidal variations at Onsala are relatively small, reaching only 20-30 cm in amplitude. However, being in a quite shallow bay, the water level is sensitive to changes in atmospheric pressure and thus the sea surface can experience overall tides that are approximately 1.5 m peak to peak (Löfgren et al. 2014).
The second installation, SPBY in Spring Bay, Australia, differs from GTGU in that it has never been intended as a GNSS-R installation. However, because of its location on a pier, it happens to have a clear view of open water and can thus be used for GNSS reflectometry. The view is, however, narrower than the one of GTGU because of other piers in the area, as presented in the bottom panel of Fig. 4. Again, the installation uses standard geodetic equipment with a LEICA AT504 choke ring antenna and a LEICA GRX1200, also recording GPS and GLONASS to RINEX files with 1 Hz sample rate.
Colocated with the GNSS installation, there is also a tide gauge measuring the sea level using an acoustic Aquatrak sensor and a Vegapuls radar sensor. The data are recorded at 1 Hz but reported as 1-min averages. Lastly, the tides are also different in characteristic in Spring Bay, with amplitudes ranging up to 0.9 m.
As described by Strandberg et al. (2016), the data collected from the two GNSS installations are discriminated depending on elevation and azimuth using a manually constructed, site-dependent mask so that only reflections from the area of interest are considered in the processing. No discrimination is made between GPS and GLONASS satellites, and both L1 and L2 frequencies are used.

Measurements of sea surface height with Kalman filtering
One feature of the combination of B-splines and the Kalman filter that is of particular interest is that the solution can be evaluated in both real time and with a delay. This is because the nodes i−3 , … , i affecting the measurement at a time t continue to be updated until i has been shifted out of the state vector. Therefore, there can be some delay until the final solution is presented, depending on the temporal spacing of the B-spline nodes. However, the sea level can of course be evaluated at any point during this time, even directly after the measurement at time t .
Here, we present the performance for both real time and final solution at our two test installations and compare it to other GNSS-R techniques. The final solution is, of course, interesting for historical records and long-term monitoring of the sea level and its change derived from time series, see for example Fig. 5, whereas real-time solutions are of interest for applications in need of direct measurements, for example marine shipping (Pugh 2004), reservoir monitoring (Becker and Yeh 1974), sea state forecasting and warning systems (Holgate et al. 2008), and similar applications. Of the two techniques, we compare with here, spectral analysis with height-rate correction and least-squares inversion, only the former is close to being in real time. Therefore, the real-time performance is primarily compared to the performance of that method, whereas the final solution will also be compared to the more precise inverse modeling method.
For both SPBY and GTGU, the comparison statistics are calculated from an 80-day period of measurements, from day 5 to 85 of the year 2018 for SPBY and from day 20 to 100 of the year 2017 for GTGU. Subsets of the tide gauge data and GNSS-R solutions from the respective site are presented in Fig. 5. During the data collection period, the maximum tidal difference measured by the tide gauges at both sites was similar, at 1.52 m for SPBY and 1.36 m for GTGU. However, as is clearly visible from the figure, the dynamics are different at the sites, with sub-daily tides dominating the sea level at SPBY, whereas other more irregular effects are driving the sea-level variations at GTGU. Therefore, the rate of change at SPBY is generally higher, degrading the performance of the spectral methods in particular because of the underlying assumption of static surfaces and the limits of the first-order correction. In addition, as mentioned in the previous section, the area from which data are collected at SPBY is smaller than at GTGU, which can affect the precision of the inverse modeling and Kalman algorithm, as more data and better temporal coverage have been shown to improve the solutions based on the inverse model (Strandberg et al. 2016). Figure 6 depicts the performance of the Kalman filter when evaluated at different delay times, i.e., evaluating the sea height at time t using the coefficients from the state vector x t+Δt . In order to compare the performance of the different inversion approaches, we compute the standard deviation of the retrieved GNSS-R sea-level values with respect to the colocated tide gauges which gives a measure of the precision of the measurements while ignoring systematic offsets arising from different reference levels.
The first feature to notice in Fig. 6 is that the Kalman filter is the only technique that provides results in real time. The spectral analysis method can only deliver results once a satellite has passed a large enough elevation window, which takes roughly 20 min. In general, the least-squares method has significantly longer delay times (Strandberg et al. 2016). In addition, for both SPBY and GTGU, the precision of the real-time Kalman retrievals is significantly better than for spectral analysis by a factor of four and two, respectively. This makes real-time sea height retrieval with GNSS-R much more viable for fields such as marine transportation where measurements need to be both precise enough and provided without delays.
As can be expected, the standard deviation decreases further with time as the B-spline amplitudes in the state vector are updated with more measurements. At GTGU, the iteration process leads to a final RMS error of 1.48 cm, even outperforming the inverse modeling, with an RMS error of 1.73 cm. This level of precision is actually reached already within 30 min. At SPBY, it takes slightly longer for the convergence, finally reaching an RMS error of 3.25 cm, within 10% of the least-squares precision.

Retrieving properties of the reflecting surface
As with the least-squares inversion, the height parameter is not the only parameter of interest. The amplitude, phase, and damping parameters are all affected by various properties of the surroundings. In an earlier paper by Strandberg et al. (2017), it was shown that these effects are susceptible to changes in the antenna environment and that especially the damping is useful to determine the presence of ice on the water surface and that least-squares inversion is capable of retrieving the information. In Fig. 7, we present that this is possible with Kalman filtering as well, which is to be expected as both use the same underlying model of the SNR. Like for the previous algorithm, the damping values retrieved by Kalman filtering are seen to decrease during the period when ice is reported outside of the Onsala peninsula by official ice maps (SMHI 2017).
As the Kalman filter reports instantaneous values, whereas the least-squares inversion reports daily average values, the time series presented in Fig. 7 differ slightly. The largest difference is around the transition periods, where the Kalman time series responds on a time scale of a few hours, whereas the inverse modeling has a trend over a few days. This can be seen in both the freezing and thawing at the start and end of the ice period of 2012, but also in the shorter thawing events near the end of the period, which are better resolved using the Kalman filter. With the increased time resolution of the damping retrieval, we have reached a point where now the validation data set does no longer provide the temporal and spatial resolution in order to catch the small-scale changes of ice coverage on the coast detected by GNSS-R. Upper graph: time series of data from the Onsala tide gauge and GNSS-R retrievals using data from GTGU. Lower graph: tide gauge data and GNSS-R retrievals from Spring Bay. The mean of each series has been removed to aid comparison 1 3 61 Page 10 of 12

Discussion
The performance of Kalman filters is sensitive to the process models used in the update step. For this study, the process noise settings have been selected through empirical testing. However, as the number of parameters in the model is quite large, i.e., ten or more depending on the number of GNSS systems and frequencies, the number of choices for the process noise levels of the individual parameters becomes very large, which implies that not every combination can be explored. Therefore, it is likely that the performance of the algorithm could be improved even further by adjusting the noise models. Instead of doing this through trial and error, sufficient knowledge about the underlying dynamics of the parameters involved could help to set the process noise values to a more physically correct level. However, this would require more detailed knowledge about all factors that can affect the SNR and how these vary over time. Thus, the problem of finding correct process noise level is site-specific and includes, for example, the exact geometry around the receiver, atmospheric conditions, vegetation growth cycles, and more (Nievinski and Larson 2014a). This promises to be an interesting topic for further studies.

Conclusions
We show that real-time GNSS reflectometry is possible using Kalman filtering techniques, with the unscented Kalman filter and a dynamic B-spline approach. We apply it to deriving sea surface heights, and the filter solution agrees well with measurements from colocated tide gauges, Fig. 6 Standard deviation of the difference between the tide gauges and the Kalman filter solution as a function of time before evaluation for both GTGU and SPBY. Spectral analysis is possible after recording data from a sufficiently long satellite arc, i.e., roughly 20 min (Williams and Nievinski 2017). Least-squares inversion (LSQ) requires significantly more data, and thus, time before analysis is possible. Still, for comparison, the precision of the technique is indicated in the figure using the shaded areas outperforming the previous spectral methods in real time and producing measurements on a similar level of precision as least-squares inversion in post-processing. Therefore, with this new algorithm, real-time GNSS-R can be a feasible alternative to traditional tide gauges, for example, in areas where remote sensing techniques are preferable.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.