An Evaluation of Sun-Glint Correction Methods for UAV-Derived Secchi Depth Estimations in Inland Water Bodies

1. Introduction

Secchi depth (SD), an essential measure of water transparency in aquatic ecosystems, provides a critical indication of water quality and ecological health [1,2,3]. In Europe, SD is important for following the water quality rules set by the European Water Framework Directive (WFD), which mandates member states to uphold ecological standards [4]. SD thus helps identify impairments affecting water quality and aids in devising and executing management plans to protect water bodies [5].

The Secchi disk is the main tool used for obtaining SD measurements, particularly in vast and complex aquatic ecosystems [6,7,8]. However, traditional methods – visual observations using a white Secchi disk, can be labour-intensive, time-consuming, and potentially influenced by observer bias [9,10]. They also typically offer point-level measurements, which may not reflect water transparency over larger areas [11,12].

The limitations of traditional methods have spurred interest in remote sensing techniques for estimating SD in complex aquatic ecosystems. Techniques ranging from satellite data to Unmanned Aerial Vehicles (UAVs) equipped with Light Detection and Ranging (LiDAR) systems and hyperspectral imagery have been explored for accuracy [13]. These techniques promise cost-effective and efficient high-resolution spatial data on water transparency over large areas, with multispectral imagery from UAVs emerging as a promising alternative [14].

SD is primarily influenced by three optical components: chlorophyll-a (Chl-a), coloured dissolved organic matter (CDOM), and total suspended matter (TSM) [15,16]. Various algorithms exist to calculate water parameters from multispectral data, including mechanistic models [17], artificial neural networks [18], and regression algorithms [19]. Each offers unique advantages and applicability, depending on the specific characteristics of the data and the aquatic ecosystem under study. For instance, a study by Chusnah and Chu [20] demonstrated the application of machine learning in estimating chlorophyll-a concentrations, which are commonly used as indicators for assessing the trophic level of lakes and the state of water quality. The study utilized machine learning to implement a band ratio algorithm and generate chlorophyll-a maps from Sentinel-2 and Sentinel-3 satellite images. However, mechanistic models, relying on physical laws and principles to simulate light attenuation and scattering processes in water, often provide the most reliable and accurate predictions, particularly when dealing with smaller datasets, that would not be sufficient enough for machine learning models [21].

Lee et al. [22] provided a foundational understanding of the optical properties influencing SD, which significantly contributed to the development of the Quasi-Analytical Algorithm (QAA) for more accurate and reliable water clarity estimations. The QAA, a commonly used mechanistic model for SD, has been utilized in various water bodies [21,23,24]. It has been applied to MODIS and MERIS satellite data, where it reduced the root mean square error (RMSE) of SD estimation from 1.5 m to 1.0 m. Furthermore, the QAA has been used to account for the residual error in reflectance data from MODIS satellite data, demonstrating its potential for remote sensing in monitoring and managing water resources [25]. The algorithm showed excellent results (R² = 0.96, MAPD = 0.18) when validated with independent measurements covering oceanic, coastal, and lake waters [26]. However, previous studies have not delved into the potential benefits of utilizing multispectral cameras onboard UAVs. The integration of these cameras with UAVs offers potential improvements in spatial resolution and data availability, bridging the gap between in situ and satellite remote sensing measurements [27].

It is important to accurately account for reflected light from the water surface, more specifically the sun-glint effect, as it can lead to inaccuracies in further processing of water quality algorithms [27,28,29]. The most simple way of avoiding sun-glint is careful UAV flight time and direction planning, however, since the water surface is often uneven and it is hard to reduce sun-glint completely [30]. There are several methods of reducing the sun-glint effect in multispectral UAV images during post-processing, for example, M Muslim et al. [29] tested four methods proposed by Lyzenga et al. [31], Joyce [32], Hedley et al. [33], and Goodman et al. [34] and applied them to either the whole image or just the glinted area. Other studies used methods of Hochberg et al. [35], removing glinted pixels as NIR threshold, or using HydroLight simulations [27]. The application of these methods should align with the specific requirements of the study. For instance, if the primary objective entails bottom mapping, the method proposed by Lyzenga offers superior results, as demonstrated in M Muslim's et al. [29] study. On the other hand, if the analysis is focused on assessing surface water quality parameters, the Hedley method emerged as the preferred choice by Windle and Silsbe [27]. Thus, the selection of sun-glint correction techniques requires careful consideration of the objectives, ensuring the most effective and accurate outcomes in different scenarios.

This research aimed to evaluate the effectiveness of QAA for the determination of SD using multispectral cameras onboard UAVs. It was done by testing QAA on image datasets pre-processed using five different methods of sun-glint correction. Additionally, the study examined how water constituents and solar zenith angle affected the discrepancy between actual and modeled SD values. We hypothesise that the radiometric accuracy of calibrated UAS sensor should meet the required accuracy of 5 %, expected with ocean colour remote sensing when compared to in situ measures of Hooker et al. [36] and application of sun-glint correction methods will improve usability in SD modeling, by deviation of RMSD of the same 5 % from in situ measures, that would still provide practical results for water quality assessment.

2. Materials and Methods

The The study was performed in 43 water bodies in Lithuania (Figure 1), represented by high variability and proportion of optically active in-water components – turbidity, Chl-a and CDOM (Table S1), thus bolstering the robustness and generalizability of the findings. The research was conducted from May 2021 to May 2023.

2.1. UAV Data Collection and Calibration

Inspire 2 Unmanned Aerial Vehicle (UAV) was used to facilitate the acquisition of data, equipped with a RedEdge-MX camera [37]. The camera features five bands: blue (475 nm ± 16 nm), green (560 nm ± 13 nm), red (668 nm ± 8 nm), red edge (717 nm ± 6 nm), and near-infrared (NIR) (842 nm ± 28 nm). To ensure optimal image capture, the UAV was flown at a height of 60 m to optimize the balance between spatial resolution and the area covered in each image allowing for both relatively high-resolution imagery (from 3 to 4 cm/pixel) and a reasonable area coverage per flight (around 40 – 50 m2), with the camera programmed to capture images every three seconds. This frequency was necessitated by the inability to view the live feed and capture images in desired areas from the Micasense RedEdge-MX camera at a distance. The maximum distance between the SD measurement site, where GPS coordinates were recorded, and the image GPS coordinates were approximately 40 m, while the minimum distance was 0 when the measurement was above the boat. The operation of the UAV was coordinated by two individuals — one controlling the UAV from the shore and the other onboard a boat, both communicating via phone during the Secchi disk measurements. The onboard operator also concurrently measured other water quality parameters outlined in Section 2.2.

For image calibration, the original Micasense RedEdge calibration panel was used. According to Micasense's guidelines, a calibration panel image was captured at the shore before each flight, ensuring no shadows were cast on the reference panel or drone [37]. Three images were taken each time for assurance.

2.2. In Situ Data

Concurrent with the UAV image capturing, in situ SD measurements were performed using a 30-cm white Secchi disk, a widely accepted method for assessing water transparency in lakes, rivers, and oceans [38].

The parameters of Coloured Dissolved Organic Matter (CDOM) – the light-absorbing component of dissolved organic matter, Chlorophyll-a (Chl-a) – an indicator of phytoplankton biomass, and turbidity – a measure of the cloudiness or haziness of water were measured alongside SD to assess their potential influence on water transparency.

Water samples for Chl-a measurements were filtered through glass fiber GF/F filters with a nominal pore size of 0.7 μm and extracted into 90% acetone. Photosynthetic pigments were measured spectrophotometrically and estimated according to the trichromatic method [39,40]. CDOM was measured spectrophotometrically after filtration through 0.22 μm membrane filters. The CDOM absorption coefficient at 440 nm was derived according to Kirk [41]. SHIMADZU UV-2600 spectrophotometer was used for the analysis of Chl-a and CDOM. Turbidity was measured with a portable turbidity meter (Eutech InstrumentsTN-100, Landsmeer, The Netherlands) in the Nephelometric Turbidity Unit (NTU). The instrument has a light-emitting diode in the near-infrared range (Hach Lange at 860 nm and Eutech Instruments at 850 nm), and the detector measures the scatter at a 90° angle. This method is based on the International Standard Organisation (ISO) 7027.

The average SD of the lakes was 1.91 ± 1.53 m (± – standard deviation), ranging from 0.25 to 7.2 m. The mean CDOM was 2.95 ± 3.2 m⁻¹ (min-max: 0.37 – 20.01 m⁻¹), the mean Chl-a concentration – 26.59 ± 26.48 (1.13 – 113.23) and the mean turbidity – 6.76 ± 11.8 NTU (0.00 – 70.62 NTU).

In situ remote sensing reflectance Rrs was acquired to validate wavelengths from UAV observations of water surface reflectance. Rrs was measured in the spectral range of 400–800 nm by the simultaneous measurements of downwelling irradiance, upwelling radiance, and downwelling radiance, performed with a WISP-3 spectroradiometer [42]. Rrs was calculated according to equation 1:

$$R_{rs}=\frac{L_u-\rho L_d}{E_d}$$

(1)

where Lu is the upwelling radiance, Ld – the downwelling radiance, Ed – the downwelling irradiance, ρ – a water surface reflectance factor equal to 0.028.

Central wavelengths of the Rededge MX camera were used (475 nm, 560 nm, 668 nm, 717 nm, 842 nm). NIR was excluded from comparison with in situ reflectances, which were not measured at 842 nm, however, NIR was still used for the sun-glint correction step section 1.4.

2.3. Pre-Processing the UAV Data

The collected UAV images underwent pre-processing of correcting atmospheric effects and standardizing the format suitable for analysis (Figure 2). The conversion of raw images to radiance involves several corrective steps to eliminate biases and errors, accounting for dark pixel offset, vignette effect, as well as aligning images due to distances between sensors, which could potentially affect the accuracy and reliability of the data. These steps were achieved using the Micasense Python (version 3.7, Python Software Foundation, 2018) workflow, mainly the function of raw_image_to_radiance, as described in the Micasense Github repository for users [37].

Since reflectance is relatively independent of illumination conditions [11,43], it is preferred for calculating SD using remote sensing images. The conversion of radiance to reflectance was performed using a reference panel of known reflectance to determine a scale factor between radiance and reflectance. This scale factor was applied to the entire image to obtain a reflectance image. The accuracy of the reflectance image was verified by extracting and checking the same reflectance panel region used to calculate the scale factor for any trends or inconsistencies [37].

Reflectance was further normalized by dividing it by π, assuming water as a Lambertian body, following the mathematical model (2).

Reflected radiance = ρ * (cos θ_i)/π

(2)

Here, cos θ_i represents the cosine of the angle between the incident light and the surface normal and ρ – a surface reflectance factor [43]. This normalization ensures that the reflectance values fall within a standardized range (0 to 1), facilitating consistent comparisons and calculations across different surfaces, lighting conditions, and measurement devices [44,45].

2.4. Image Masking and Sun-Glint Correction

From the initial set of water bodies, four were eliminated (Tūbausių Reservoir, Grabuostas Lake, Ilgis Lake, and Mušėjus Lake) from further analysis. One was excluded due to excessive cloud glint, while the other three measurements were performed too close to the shore or contained emerged macrophytes, leaving no sufficient area where the reflectance would not be affected (Figure 5a). One water body (Kruminių reservoir) was surveyed at a late time when the zenith angle was approximately 80 degrees, resulting in shadows covering a half portion of the area; however, the unshaded area could be used for the analyses (Figure 5b). In total, 39 water bodies were left after discarding unsuitable ones for statistical analysis.

Figure 3. (a) RGB of an image affected by shadows, bottom reflectance and floating macrophytes after thresholding (Grabuostas Lake), (b) RGB of contaminated image (Kruminių reservoir) with shadow and shore but where values can still be used by selecting an area (red square) that was good, (c) RGB example of a relatively good image (Musia Lake) before removing sun-glint and boat.

Figure 3. (a) RGB of an image affected by shadows, bottom reflectance and floating macrophytes after thresholding (Grabuostas Lake), (b) RGB of contaminated image (Kruminių reservoir) with shadow and shore but where values can still be used by selecting an area (red square) that was good, (c) RGB example of a relatively good image (Musia Lake) before removing sun-glint and boat.

To ensure accurate SD measurements, objects (boats, coast, macrophytes) and shadows potentially affecting the measurements' accuracy were removed. Firstly, the Normalized Difference Water Index (NDWI) was calculated to distinguish between water and non-water pixels (3). This index was then classified into binary parts using Yen thresholding [46,47], where non-water pixels were masked as NaN values. Yen method was used also on just the NIR band to classify areas where pixels were affected by sun-glint and masked as NaN in one of the sun-glint removal tests. This sun-glint removal method was later called the threshold-removed glint method.

$$NDWI=\frac{RrsGreen-RrsNir}{RrsGreen+RrsNir}$$

(3)

In 11 images water waves caused a substantial amount of sun-glint (Figure 4a), despite images being captured early in the morning around 10 am GMT+3 when the sun reflection from calm water should not have reached the lens. Therefore, the sun-glint correction was performed using the algorithms developed by Hedley et al. [33], Goodman et al. [34], Lyzenga et al. [31], and Joyce [32], which assume that the water surface reflectance is a linear combination of water reflectance and sun-glint reflectance. The models were fitted using a set of training data consisting of image pixels where sun-glint was absent, calculated from an area with 10 % lowest value NIR pixels [27]. The model was then used to predict water reflectance for sun-glint-affected pixels, reducing reflectance values according to equations (4, 5, 6, 7).

The Hedley method (4) calculates the remote sensing reflectance (Rrs) for each pixel in each band. The method first finds the minimum NIR value and then calculates the slope for each band using linear regression. The Lyzenga and Joyce methods (5 and 6) are similar to the Hedley method but use respectively the mean and mode of the lowest 10% of NIR values instead of the minimum NIR value. The Goodman method (7) calculates the Rrs for each pixel in each band using a constant A (0.000019) and a factor B (0.1) that was multiplied by the difference between the red and NIR bands.

Hedley Rrs_λ = λ + bi_λ ∗ (NIR - NIRmin)

(4)

Lyzenga Rrs_λ = λ + bij_λ ∗ (NIR - NIRmean)

(5)

Joyce Rrs_λ = λ + bi_λ ∗ (NIR - NIRmode)

(6)

Goodman Rrs_λ = λ – NIR + δ, where δ = A + B ∗ (Redλ - NIRλ)

(7)

where λ is the band of interest (blue, green, red, red edge, NIR), bi_λ is the slope of the band, bij_λ is equal to covariance between λ and NIR divided by variance in NIR, and NIR_min NIR_mode NIR_mean is the minimum, mode or mean NIR value respectively. For Goodman’s method, δ is a constant offset across all wavelengths where A and B are constants (A = 0.000019 and B = 0.1). Additionally, removing only sun-glint affected areas (determined by the binary thresholding Yen algorithm) was compared as an alternative approach to sun-glint correction methods (Figure 4c).

The images of four water bodies (Vembutų, Stebuliškių, Pakapės, and Krūminių reservoirs) were also affected by either shore or cloud glint artefacts. However, most of these artefacts were successfully removed using binary thresholding, except for one image (Figure 3b) where the area had to be manually chosen because the thresholding method was unable to accurately separate the unwanted areas.

2.5. Secchi Depth Model

Quasi-analytical algorithm (QAA) proposed by Lee et al. [22], provides a robust framework for monitoring SD (8), particularly in scenarios where in situ measurements may be unavailable. Compared to the empirical approach, this semi-analytical method offers a significant advantage – it does not necessitate the recalibration of the retrieval algorithm with in situ data [23]. This enhances its utility in diverse settings. As such, it was considered more suitable for monitoring SD in various water bodies.

Lee et al. [22] paper introduced a mechanistic model that accounts for the effects of light attenuation, scattering, and reflection in the water column, as well as the properties of the Secchi disk itself, in determining the SD. The parameters that determine SD in this algorithm are the total absorption coefficient 'a and the total backscattering coefficient 'bb. From these parameters, the diffuse attenuation coefficient Kd was calculated, which is the main variable in the SD formula, besides the Rrs band.

$${\mathrm{K}}_{\mathrm{T}}\text{\_}{\mathrm{K}}_{\mathrm{d}}=\frac{1.04\ast {\left(1+5.4u\right)}^{0.5}}{1/{\left(1-\frac{sin{\left(\theta \right)}^2}{R{I}^2}\right)}^{0.5}}$$

(8)

$$\mathrm{SD}=\frac{1}{kt\_kd \ast minKd} ln\left(\frac{0.14-minRrs}{0.013}\right)$$

(9)

where minKd – minimum value Kd chosen out of Kd calculated with blue, green and red band, minRrs – above surface remote sensing reflectance of band that had the lowest Kd value. Kt_Kd – upwelling radiance diffuse attenuation coefficient and was used instead of a constant value of 1.5 as suggested by Jiang et al. [23], where θ is the solar zenith angle, RI is the refractive index value of water equal to 1.33 and u denotes the ratio of backscattering coefficient to the sum of absorption and backscattering coefficient.

The sun zenith angle was computed from in situ measurements using the time and location of the observation. This information was then used to calculate the solar position using the Python library pytz [48]. This approach provided a reliable means of determining the solar zenith angle, which is an important parameter in applications related to SD acquisition from Rrs.

Selecting the appropriate reference wavelength is of major importance for the final SD value [26]. The green band (560 nm) was used as a reference if the red band Rrs (668 nm) is < 0.0015 sr⁻¹, otherwise the red band (668 nm) was used, with accordingly modified calculations of other parameters, as suggested by Lee et. Al. [22].

SD was calculated using the QAA method for images of 6 types: corrected after Hedley, Lyzenga, Joyce and Goodman sun-glint algorithms, threshold-removed images and images with no correction.

2.6. Validation and Interpret the Results

All of the final modeled SD values and reflectance images were averaged to ensure that images can be compared with in situ data. The final reflectances obtained after applying various correction methods and no correction (1.4 – 1.5 sections) were validated with in situ reflectances (WISP-3). The accuracy of these methods was evaluated by bias (9), while the root mean square deviation (RMSD) was used as an indicator of the QAA model‘s precision (10), and the Pearson’s correlation coefficient (r) described the relationship strength between the model‘s output and the real data values.

$$\mathrm{Bias}=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N\left(X_{estimated,i}-X_{measured,i}\right)$$

(9)

$$\mathrm{RMSD}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(X_{estimated,i}-X_{measured,i}\right)}^2}{N}}$$

(10)

The same accuracy and precision measures were applied to modeled SD values. The Generalized additive models (GAMs) were employed [49] to investigate the relationships between the difference in modeled and in situ SD (the response variable) and a set of independent variables: CDOM, Chl-a, turbidity, and Solar Zenith. The GAMs were chosen due to their flexibility in modeling non-linear relationships and their ability to handle interactions between predictors. Before analysis, the crosscorrelation (based on the Pearson correlation coefficient) between the independent variables was determined. The correlation was relatively high (r = 0.76) between Chl-a and turbidity, therefore turbidity was not included in the GAMs. Before the interpretation of the GAMs results, the residuals were visually inspected with diagnostic plots for normality and equal variance against the fitted values. F and p values were obtained to assess the relative importance and significance of the independent variables (p < 0.05 was considered as a statistically significant relationship). The fit of the model was evaluated using the explained deviance. A response plot was graphically represented to visualize the relationship between an explanatory variable and the response. As there were multiple predictors, each one was plotted separately with a smooth curve and the confidence interval of the effect.

To visualize the SD results, a subset of images with calculated SD values were exported as TIFF and mosaiced according to image metadata GPS coordinates in QGIS version 3.16 [50]. The result was then visualized on top of the RGB mosaic that was mosaiced using OpenDroneMap [51] photogrammetry software.

3. Results

3.1. Band Validation after Sun-Glint Correction

Across all correction methods and the in situ data, the green band (560 nm) consistently showed the peak mean value, while the blue band (475 nm) indicated the valley or lowest mean value. The general shape of the data appears to peak at the green band, with decreasing values on either side at the blue and red bands (668 nm), and a slight increase at the red edge band (717 nm). This pattern was consistent across all the correction methods and the in situ data, indicating the robustness of this spectral feature in the multispectral UAV image data.

The highest mean values across all bands were observed with the Lyzenga correction method, while the lowest was observed with the Goodman correction method (Figure 6). The in situ values were generally lower than the corrected values, with the exception of the Goodman correction method.

Figure 5. Spectral scatterplots demonstrating the agreement between in situ measured and sun-glint corrected (or not corrected) reflectances in the UAV images. Statistics (r, RMSD, bias) and points are colour coded accordingly except for the red edge that is grey. A diagonal black line marks perfect agreement (1:1).

Figure 5. Spectral scatterplots demonstrating the agreement between in situ measured and sun-glint corrected (or not corrected) reflectances in the UAV images. Statistics (r, RMSD, bias) and points are colour coded accordingly except for the red edge that is grey. A diagonal black line marks perfect agreement (1:1).

The validation results provided robust evidence regarding the performance of the correction methods and the reliability of the final reflectance values. The correlation for all methods (Figure 7) had the same trend where it was the lowest for the blue band (r = 0.38 ± 0.12) as well as for the red band (r = 0.56 ± 0.14) but higher for the green band (r = 0.75 ± 0.11) and red edge band (r = 0.80 ± 0.10). However, RMSD was relatively low for the blue band (RMSD = 0.0043 ± 0.0016) and red (RMSD = 0.0048 ± 0.0018) bands and slightly higher for the green (RMSD = 0.0060 ± 0.0013) and red edge (RMSD = 0.0053 ± 0.0014) bands it was relatively high. Bias for most methods followed a similar trend of larger underestimation for the blue band and green bands, then slight overestimation for the red band with Goodman, Hedley and threshold-removed glint methods but still undervaluation for other methods, and slightly lower undervaluation for the red edge band. Goodman‘s method was exceptional to these trends as green, red and red edge bands were overestimated.

Overall Goodman’s algorithm showed the highest correlation for the green and red edge bands, with r values of 0.90 and 0.92, respectively. However, it had low r values for the blue and red bands, with 0.32 and 0.71, respectively. The RMSD values for Goodman were relatively low across all bands, ranging from 0.002 to 0.005. Similarly, the Hedley algorithm consistently performed well across all bands, with r values ranging from 0.58 to 0.88. The RMSD values were similar to those of Goodman, ranging from 0.003 to 0.005.

The Joyce and Lyzenga algorithms showed matching performances, with lower r values and higher RMSD and bias values compared to Goodman and Hedley. The r values ranged from 0.27 to 0.69 for Joyce and from 0.26 to 0.67 for Lyzenga. The RMSD values ranged from 0.006 to 0.007 for both algorithms.

The control group with no correction applied showed moderate r values, ranging from 0.41 to 0.81. The RMSD values were similar to those of Goodman and Hedley, ranging from 0.005 to 0.006.

3.2. Validation of QAA SD Model

The performance of the QAA model‘s (Figure 8) ability to predict SD when compared with in situ measurements, showed a relatively high correlation across all methods with r ranging from 0.74 (threshold-removed glint) to 0.92 (Hedley glint correction). RMSD (from 0.65 to 1.05 m) and bias (from -0.78 to 0.58 m) showed acceptable results over all methods as well.

Comparing methods between themselves, the accuracy of all parameters had similar trends as accuracy for band comparison with in situ measurements (section 2.1), where the SD values were overestimated for the smaller SD values and underestimated for the larger ones, except for the Goodman’s method where the most of the SD values were overestimated.

Hedley’s sun-glint corrected images achieved the best results according to the RMSD measures (0.65 m), also r was just slightly smaller (0.91) than with the Goodman’s method (0.92) that achieved the best results based on the r value of 0.92. However, the RMSD value for Goodman’s method was relatively high (1.00 m) to the values of other methods. This method overestimated most of the values both for small and large SD values.

The Joyce and Lyzenga methods showed similar results between themselves with r values of 0.87 and 0.85, respectively. The RMSD values for Joyce and Lyzenga were 0.79 and 0.86 m, respectively and biases showed underestimation.

The worst performing method was when only sun-glint-affected pixels were removed and all other pixels were left unchanged (r = 0.75, RMDS = 1.05 m, bias = 0.13 m). The control group with no correction showed slightly better results: r value of 0.89, the RMSD value of 0.74 m and underestimation according to bias.

3.3. Relation with Water Constituents

The environmental factors in GAMs significantly explained (38.3 %) the variance of the difference in SD measurements. The interaction term of the sun zenith angle and CDOM was significant (F = 6.808, p-value < 0.05), suggesting that these factors together affect the accuracy of SD retrieval. The most important and statistically significant factor was CDOM (F = 10.47, p < 0.05), followed by the solar zenith (F = 4.84, p = 0.02). Chl-a values did not have a significant effect for the GAM model (F = 0.295, p = 0.59).

The modeled SD values were overestimated (to over 1.3 m) when CDOM values < 7.5 m⁻¹ (Fig. 12), while underestimated (over 1 m) when CDOM was > 7.5 (although underestimation did not significantly change when CDOM >12). For the solar zenith, the SD values were overestimated (over 0.7 m) at the lower solar zenith angles (< 45 degrees), and then a relatively low effect (within ± 0.5 m) was between 45 and 75 degrees of the solar zenith angle and the underestimation sharply increased to over 1.5 m when the solar zenith angle was > 75 degrees.

4. Discussion

4.1. Advancements in SD Measurements

This research focused on evaluating the effectiveness of the QAA for determining SD using multispectral camera onboard UAVs. This novel approach expands the current understanding of measuring SD by introducing the potential use of UAVs, balancing the broad coverage provided by remote sensing methods with the high accuracy characteristic of in situ measurements, therefore serving as a more comprehensive alternative to traditional methods.

The experimental results did not uniformly support our hypothesis that the application of sun-glint correction methods would enhance the utility of multispectral images in SD modeling by reducing the RMSD by about 5 % and reducing the bias by a similar amount. The precision observed across all bands appears to correspond with the spectral band reflectance intensity, where bands demonstrating relatively lower reflectance intensity (e.g. blue band) typically displayed a lower compliance between the reflectance measured from the UAV and the in situ reflectance measurements and the reverse held true for bands with relatively higher reflectance intensity. Environmental factors like weak water surface signal and roughness of the water surface can introduce systematic and random errors, respectively, into water surface detection [13]. The correlation with the blue were identified as the least robust (Figure 7), which can be ascribed to their increased vulnerability to scattering, a trait inherent to these bands in the water environment [52].

Among the tested glint correction methods, the Hedley and Goodman methods emerged as the most effective across all bands, with Hedley demonstrating the lowest RMSD across all bands. For instance, a study conducted by M Muslim et al. [29] employed a similar methodology to our study, testing multiple sun-glint correction methods. Their findings indicated that the Lyzenga method yielded the most accurate results. However, it's important to note that the primary objective of their research was to map coral reefs, and in most of the study area, the bottom was visible, which may have influenced their results. In contrast, Windle and Silsbe [27] found that the Hedley method provided the highest correlation coefficient (r) and the lowest RMSD, corresponding closely with the results of our study.

The key distinction between the three glint correction methods (excluding the Goodman method) lies in the slope or covariance index of NIR band values. In the case of the Hedley method, the slope was the largest, as it takes the minimum value of the band as the starting point, as opposed to the mode (Joyce method) or the mean (Lyzenga method) [31,32,33]. This suggests that the images used in our study required a larger numeric correction due to the lowest initial reflectance values.

The results showed that the Goodman's method reflectance intensity values tend to be significantly different from other methods for the blue, red and red edge bands. This method, which was ranked as the second best performing in our study, adopted a different approach, using constant values to correct for the sun-glint effect [34]. This approach led to overall better r values between in situ and the green, red, and red edge bands on-board UAV. However, the performance of the blue band was significantly lower than when using the Hedley method. Despite this, the Goodman method demonstrated the lowest RMSD values and small overestimation bias, underscoring its effectiveness in applications where low reflectance intensity deviation from in situ radiometer is required.

The results also indicated that applying glint correction to the entire image rather than just the binary thresholded area can yield better results in terms of the overall accuracy of the reflectance values. This observation implies that the method of applying glint correction solely to the binary thresholded area might not be the most effective strategy, as it may overlook potential glint effects present beyond this area. Additional sun-glint correction research should be considered in the future for better generalization of reflectance correction [53].

The overall agreement between the in situ SD measurements and the modeled SD values might be connected to their handling of in-water constituents such as Coloured Dissolved Organic Matter (CDOM) and solar zenith angles, which significantly influenced the accuracy of the models. CDOM was found to predominantly influence the discrepancies between in situ measurements and modeled SD. Given that CDOM primarily absorbs light in the ultraviolet (UV) and blue regions of the spectrum, resulting in relatively low acquired signal by sensor [54], this agrees with our observation that the blue band demonstrated one of the weakest correlations between in situ and UAV-derived reflectance measurements [55]. It is plausible that Goodman's and Hedley’s methods more effectively managed this factor where Hedley’s had the highest r and Goodman’s had the lowest RMSD for blue band reflectance, resulting in a closer match with in situ values.

Similarly, high solar zenith angles above 70 ^o impacted the accuracy of modeled SD, leading to underestimations by up to 1.5 m. This is likely due to increased scattering and absorption of light at higher zenith angles, resulting in less light reaching the water's surface and thus larger differences between modeled and in situ SD values [56]. This reinforces the need for sophisticated algorithms that can accurately model these complex environmental factors in SD predictions.

4.2. Practical Applications

Accurate SD determination in large areas is particularly important given the increasing demand for high-resolution data on water transparency for applications including water resource management, environmental monitoring, and ecological modeling [57]. Measures of SD can reduce the need for boat measurements in lakes, also allowing for data collection at a higher frequency, surpassing traditional monthly monitoring intervals that may be inadequate for dynamic water bodies with recurring algal blooms, which in some cases can be inaccessible due to terrain or vegetation around the lake [58]. As well, UAVs could be used in shallow coastal waters (< 1.5 m), where a research vessel usually cannot access them, e.g. in the Curonian Lagoon [59].

The optimal fit model using Hedley sun-glint correction on the whole image was used to construct a mosaic from 45 individual images (Figure 9) of Kašučių Lake 2021 09 20. All of the pixels were left unmasked in order to show how shadows affect the final results, and therefore this aspect should be taken into consideration while planning the acquisitions.

This example shows the efficacy of the QAA best-fit model in transforming discrete data into coherent, large-scale representation of SD. The SD varied from 4 to 6 m, where the highest values were determined in the areas covered by tree shadows. Elsewhere SD values were more homogeneous (4.0 to 4.8 m), but some noise was still apparent in the orthomosaic (snow-like effect) as SD values are calculated for each pixel. For this reason, the mean value of the whole image was compared to in situ measurements in the area, instead of selecting a small square or point in the images, which would have resulted in lower accuracies.

Overall, UAVs can provide higher-resolution images than satellites, allowing for more precise measurements of water bodies, including small and shallow ponds [60,61], which are often not well captured by satellite imagery [61,62]. Furthermore, UAVs can be flown at a specific time and location, which is particularly important for SD measurements, as the measurement of water transparency can be influenced by various environmental factors, such as partial cloud cover, water turbidity, and the angle of the sun [27,63]. Although weather conditions influence the ability of UAV flight (with Inspire 2 UAV if winds exceed around 10 m/s), they can collect data during fully cloudy conditions, where optical satellite data are heavily affected by clouds. Additionally, they can be equipped with diverse sensors for capturing multi-spectral or hyperspectral images, thus providing valuable data on water quality and ecology [14,64]. Furthermore, UAVs can serve as ground truth to validate satellite data [65,66], thereby enhancing its accuracy and monitored area, e.g. one tile coverage of Sentinel-2 can reach from 100 000 ha, while UAVs cannot cover relatively large areas.

This method has potential applications in shaping policies and regulations related to water bodies. By yielding accurate, timely, and spatially inclusive data on water transparency, it equips decision-makers with invaluable resources, necessitating immediate intervention or remediation measures. Moreover, UAV technology presents a strategic advantage in regions where conventional methodologies struggle due to logistical complications. This is particularly evident when we consider the capabilities of UAVs for mapping vast areas. A single UAV flight, which takes approximately 25 min., can effectively map an area as large as 25 ha. To put this into context, Lithuania has around 6000 lakes, and about 340 of them are larger than 5 ha [67], which would be suitable for monitoring using Sentinel-2 satellite, assuming that the shape is not elongated. The rest of the smaller lakes could be monitored by applying methodology from this study and potentially could improve the accuracy and coverage of current monitoring, which covers around 80 lakes every 6 years in Lithuania. With the methods provided in this paper, processing workflow for one 5 ha lake should not take longer than 20 min. Moreover, this process can be fully automated, requiring only the supervision of the final results.

4.3. Future Research and Potential Limitations

The QAA model possesses the flexibility to process remote sensing data sourced from an array of sensors, including but not limited to MODIS, MERIS, OLCI, MSI and GOCI. By integrating a more diverse set of global in situ measurements corresponding to various water types, it's conceivable to further refine the precision and effectiveness of the QAA model.

While the research results demonstrate potential, it is essential to emphasize that the investigation was conducted across a broad range of unique aquatic environments and included measurements captured under variable sun zenith angles. Consequently, for future work, there's potential to modify the existing QAA used for SD estimations to better account for CDOM and sun zenith angle, that in this study have shown to be critical parameters. This modification could improve our understanding and predictions in the context of diverse and changeable aquatic conditions. There are several limitations that need to be addressed in future research. One of the main limitations is the assumption of linearity in the Hedley method. This method requires at least one dark pixel unaffected by glint [33] to be visible in the image, and if the whole image is affected by glint, this deglinting method will provide incorrect data. Understanding these limitations aids in contextualizing the results and fosters the development of more refined glint correction methodologies in the future. This is also relevant for a partial cloud glint, when clouds are reflected from the water surface in some, but not all areas of the image, for the analysis of multispectral drone imagery [33,68]. Most methodological studies for water surface mapping using UAVs suggest careful planning of flight time, preferably during conditions of clear skies and low sun-glint [30], however, this restricts one of the main advantages of UAVs – getting data on demand, especially if there is a need to visit several inland water bodies at the same day.

Many studies underscore the importance of choosing the appropriate quasi-analytical algorithm (QAA) for different optical water types [13,21,26]. The optical properties of coastal and inland waters, which are primarily influenced by the concentrations of suspended particulate matter, phytoplankton, and dissolved organic material, exhibit spatial and temporal variability, leading to diverse optical water types [69,70]. Applying a universal algorithm in these optically complex waters often results in significant uncertainties [69,70]. This insight highlights the importance of optical water classification to enhance retrieval accuracy, as indicated by numerous studies that developed class-specific algorithms for bio-optical parameters and achieved improvements by applying optical classification in the retrieval of these parameters [70,71,72,73]. In our study, relatively high concentrations of CDOM significantly affected the optical properties of the water and thereby introducing uncertainties in the results. One possible approach to improve the accuracy of the QAA is to calibrate the algorithm for the CDOM-dominated water bodies.

In addition to the aforementioned limitations, another significant challenge to consider is the interference of high vegetation pixels in the analysis of multispectral drone imagery [74]. Vegetation and their shadows, especially when in close proximity to water bodies, can skew the reflectance measurements and thus influence the accuracy of Secchi depth estimates. Consequently, implementing strategies to exclude these pixels during image analysis can enhance the reliability of the measurements.

The applicability of water parameters extends beyond SD, allowing for the calculation of additional parameters. Prior research has demonstrated the feasibility of employing multispectral UAVs for turbidity [75,76,77], Chl-a [75,77,78,79], CDOM [77], TSS [78,80], cyanobacteria [78] and macrophytes [46,61,81]. Given the inherent scalability of the UAV-based methodology, it stands as a promising tool for extensive SD assessments and water quality surveys, thereby facilitating large-scale studies focused on water transparency. This advancement brings us a step closer to exploiting the full potential of UAV-based remote sensing for assessing and monitoring aquatic environments.

5. Conclusions

This study revealed that the accuracy of SD measurements is profoundly influenced by sun-glint correction methods employed in UAV flights. There was a consistent agreement across all methods and the in situ radiometric data, particularly for the green band, emphasizing the robustness of multispectral UAV image data. Among the tested methods, Hedley's method demonstrated superior accuracy (RMSD = 0.65 m) and precision, thereby significantly contributing to the accuracy of the UAV-derived SD data.

Moreover, findings underscored the significant role of environmental factors, particularly the CDOM and solar zenith angle, causing inaccuracies in SD measurements: a solar zenith angle > 70 ° resulted in an underestimation of up to 1.5 m in modeled SD, while CDOM > 12 m⁻¹ caused similar underestimations. Our research, therefore, supports the use of UAVs equipped with multispectral cameras as a viable method for SD determination in inland water bodies with SD of up to 7 m and lower than 12 m⁻¹ CDOM. The results point towards an approach capable of reaching a correlation as high as 0.91 and reducing the RMSD by up to 12.85 % (Hedley’s method), thereby enhancing the versatility and reliability of SD measurements.