A Semi-Empirical Anisotropy Correction Model for UAS-Based Multispectral Images of Bare Soil

: The recent developments in the performance and miniaturization of unmanned aircraft systems (UAS) and multispectral imaging sensors provide new tools for the assessment of the spatial and temporal variability of soil properties at sub-meter resolution and at relatively low costs, in comparison to traditional chemical analysis. The accuracy of multispectral data is nevertheless inﬂuenced by the anisotropic behaviour of natural surfaces, framed in the general theory of the bidirectional reﬂectance distribution function (BRDF). Accounting for BRDF effects in multispectral data is paramount before formulating any scientiﬁc interpretation. This study presents a semi-empirical spectral normalization methodology for UAS-based multispectral imaging datasets of bare soils to account for the effects of the BRDF, based on the application of an anisotropy factor (ANIF). A dataset of images from 15 ﬂights over bare soil ﬁelds in the Belgian loam belt was used to calibrate a model relating the ANIF to a wide range of illumination geometry conditions by using only two angles: relative sensor-pixel-sun zenith and relative sensor-pixel-sun azimuth. The employment of ANIF-corrected images for multispectral orthomosaic generation with photogrammetric software provided spectral maps free of anisotropic-related artefacts in most cases, as assessed by several ad hoc indexes, and was also tested on an independent validation set. Most notably, the standard deviation in the measured reﬂectance of the same georeferenced point by different pictures decreased from 0.032 to 0.023 ( p < 0.05) in the calibration dataset and from 0.037 to 0.030 in the validation dataset. The validation dataset, however, showed the presence of some systematic errors, the causes of which require further investigation.


Introduction
Diffuse reflectance soil spectroscopy provides a good alternative that may be used to enhance or replace conventional methods of soil analysis, as it is rapid, timely, less expensive, non-destructive, and allows for the simultaneous characterization of various soil properties [1]. Soil spectroscopy is currently a mature discipline thanks to continuous improvements since the mid-1960s [2]. The recent developments in unmanned aircraft systems (UAS) [3] and technological advancements in the performance and miniaturization of lightweight spectrophotometers [4] provide a toolset for the flexible assessment of the spatial and temporal variability of soil properties at a high resolution (<1 m). Hyper and multispectral imaging sensors, in particular, have been successfully employed for the UASbased remote sensing of several soil properties in several studies [5][6][7]. The accuracy of spectral remote sensing, however, is challenged by environmental conditions (illumination and atmospheric interactions), measurement protocols, data-processing procedures, and specific sensor characteristics. It is thus critical to understand both the potential and the limitations of this technology before using its data products to study earth system processes [8].
One of the biggest challenges for spectral sensing is that the reflectance of natural surfaces shows anisotropic behaviour, i.e., measured reflectance varies with illumination and viewing geometry. This results in spectral variability not only from different angular measurements and times of the day, but also within the field of view of single images. The theoretical framework for modelling the angular effects occurring between the illumination source, sensor and measured surface has been formulated as the bidirectional reflectance distribution function (BRDF) [9]. Several physical, empirical or semi-empirical models have been developed in an attempt to normalize spectral data with the parametrization of the BRDF [10]. Traditionally, BRDF model calibrations were obtained by performing multi-angular reflectance measurements using complex goniometer setups, but currently, aerial digital imaging provides an alternative platform to field goniometer systems for performing multi-angular reflectance measurements [11]. The advantage in the employment of imaging setups lies in the fact that the field of view of each image can provide a range of measurement geometries and, thus, spectral interactions between the illumination source and the measured surface [12]. In this regard, several empirical methods were developed to use the statistics of images taken on flat homogeneous surfaces to normalize off-nadir reflectance values with those obtained at nadir, with the use of an anisotropy factor (ANIF) [13,14].
Addressing BRDF-related issues in data from remote sensing and imaging spectroscopy is of primary importance to ensure their accuracy before attempting any further analysis. For a few years, UAS multispectral data have been used to infer the spatial variability of soil properties at very high spatial resolution. However, to the best of our knowledge, a thorough understanding of the bias introduced by using uncorrected image datasets is still very limited and a generally applicable workflow to construct and apply bare soil BRDF models is not available. A preliminary analysis from our own bare soil multispectral datasets that were not corrected showed that the resulting orthomosaics exhibited prominent spectral aberrances in the form of "stripes" of different reflectance intensities, which were parallel and coincided with the drone flight lines. The full potential of multispectral imaging for bare soil applications has therefore not been fully reached.
The primary objective of this article is to elaborate a semi-empirical correction methodology for multispectral imaging datasets to account for the effects of the BRDF, based on the calculation of an anisotropy factor (ANIF) calibrated from UAS-based hemisphericaldirectional [6,14] bare soil reflectance measurements. A secondary objective is to evaluate the beneficial effects of this correction method on the quality of structure from motion (SfM)-based multispectral orthomosaics, using an elaborated UAS data collection campaign of 15 flights conducted in the Belgian loam belt.
To accomplish this, sensor position and orientation data obtained from GPS measurements and SfM-based bundle adjustment algorithms were used to relate observation and illumination angles to soil anisotropy effects at the image level. The resulting data were then used to test different modelling approaches and to provide a consistent image correction workflow. Finally, the orthomosaic products obtained with commercial SfM software using original and corrected images were analytically compared for an analysis of the reduction/elimination of the "striping" artefact and the stability of the spectral measurement when explicitly considering BRDF effects in bare soils.

Multispectral Imaging System
The multispectral imaging system used to collect data was a MicaSense RedEdge-M (Micasense Inc., Seattle, WA, USA). This widely used system consists of an array of five synchronized cameras, each with a distinct narrowband filter providing spectral bands centred at 475, 560, 668, 717 and 840 nm wavelengths, with a bandwidth (full width at half maximum) of 20, 20, 10, 40 and 10 nm, respectively. The image resolution is 1280 × 960 pixels with a Remote Sens. 2022, 14, 537 3 of 15 3.75 µm pixel pitch, and the camera features a focal length of 5.4 mm, resulting in a 47.2 • cross-track field of view (FOV) and a 36.2 • along-track FOV. A diffuse irradiance sensor array with similar spectral band characteristics (the so-called downwelling light sensor, or DLS) is mounted on top of the drone in the same optical axis as the camera array and is wired to the main unit to provide synchronized irradiance data. The DLS unit also includes an on-board standalone L1 GPS and inertial measurement unit (IMU) to provide location and orientation measurements (positioning accuracy of 2-5/10 m in xy/z and orientation accuracy of 2-5/5-15 • in roll, pitch/yaw). The multispectral imaging system stores a set of five TIFF image files (corresponding to the spectral bands) for every image acquisition point, for which their respective irradiance level and a common GPS position and IMU orientation are saved in the image EXIF and XMP metadata. A calibrated reflectance panel (CRP) of ca. 50% grey reflectance provided by the camera manufacturer was used to convert raw image radiance values to reflectance.

Data Collection and Pre-Processing
A series of UAS flights with the Rededge-M rigidly mounted on a custom-made quadcopter drone was conducted over bare soil fields in the Belgian loam belt between 2018 and 2021 (flight information is summarized in Table 1, locations are shown in Figure 1a,c). The flights were conducted under clear sky conditions, on dry soil. The fields were in seedbed conditions, with very low roughness (example provided in Figure 1b), very low stoniness, no growing vegetation and very limited vegetation residue. Soils in the area are classified as silt-loam by texture according to the USDA classification. Flight altitude was ca. 45 m above ground level, resulting in ca. 3 cm ground sample distance (GSD). Flight speed was set to 5-6 m/s to ensure a frontal image overlap of 80%, while interline distance was set to obtain a lateral image overlap of at least 75%. Flight plans were designed with GIS software and then executed using the DJI GSPro application (SZ DJI Technology Co., Ltd. Shenzhen, China [15]). The camera was mounted on the UAS in a fixed position, with a frontal tilt (pitch) manually adjusted to provide an almost nadir position at the chosen flight speed and with the UAS always turning the head in the flight direction. Flight lines were oriented along the long edge of the field to provide maximum coverage efficiency, resulting in variable solar azimuth angles (given in Table 1). Table 1. General flight information of the UAS flight database: name, execution date, lateral overlap %, flight line and sunlight direction and whether the dataset was used for ANIF model calibration or validation. Names marked with an (*) represent datasets that in the original orthomosaic generation phase showed a prominent "striping" artefact effect. Flight altitude was ca. 45 m above ground level in all cases. sion to reflectance: automatically managed by PIX4Dmapper for the original orthomosaics versus manually pre-calculated on the input images for the corrected orthomosaics. In this last case, radiometric calibration, vignetting correction, black level compensation and final conversion to reflectance for the Rededge-M images were performed following the procedure described on the manufacturer's website [17] using R software (R core Team, Vienna, Austria [18]). PIX4Dmapper (PIX4Dmapper SA. Prilly, Switzerland [16]) was the SfM photogrammetric software used in this study for multispectral reflectance orthomosaic generation. PIX4Dmapper organizes the image processing workflow in three phases: (i) initial processing for key point identification and camera calibration (for both interior and exterior orientation parameters); (ii) point cloud densification; (iii) DSM (digital surface model) and multispectral orthomosaic generation. The first two phases were executed with the same processing options with both the original and the corrected image datasets. The first phase was of primary importance to obtain the adjustment of the camera orientation parameters ω, φ, κ (rotation with respect to the target projected coordinate system): since the original camera inertial navigation measurements can feature inaccuracies of up to about 5-20 • , more accurate values from a bundle adjustment algorithm were deemed necessary for the calibration of our ANIF model. The third phase differed in the choice of radiometric calibration, vignetting correction, black level compensation and final conversion to reflectance: automatically managed by PIX4Dmapper for the original orthomosaics versus manually pre-calculated on the input images for the corrected orthomosaics. In this last case, radiometric calibration, vignetting correction, black level compensation and final conversion to reflectance for the Rededge-M images were performed following the procedure described on the manufacturer's website [17] using R software (R core Team, Vienna, Austria [18]).

ANIF Correction Model
An empirical correction methodology for spectral imaging is proposed to account for the effects of the BRDF, based on the calculation of an anisotropy factor (ANIF) from drone-based hemispherical-conical [19] reflectance measurements (in the context of this study, given that the pixels of imaging spectrometers feature fairly small measuring cones, Remote Sens. 2022, 14, 537 5 of 15 their measurements are treated as an approximation of directional measurements and the resulting quantities as hemispherical-directional reflectance factors [12]). The ANIF represents the ratio between the reflectance measured at the nadir position and the reflectance measured at other angular positions in an imaging sensor. The application of this factor aims to normalize the reflectance values at each pixel position with those measured at the nadir. Since reflectance anisotropy is related to the combination of viewing geometry (sensor azimuth and zenith) and illumination geometry (solar azimuth and zenith), a linear modelling approach of the relative sensor-pixel-sun angles was tested to model the ANIF.
For model calibration, and to isolate the effects of directional reflectance on soil images, twelve flights (indicated in Table 1) were selected for the conjunction of high homogeneity of soil colour, low roughness (seed-bed conditions), and low/no slope. Within the image dataset of each flight, a subset of about 25-30 images (per spectral band) shot from the centre of each field was selected to obtain images exclusively representing flat, homogeneous, bare soil. These images served as the model calibration dataset. Radiometric calibration, vignetting correction, black level compensation and final conversion to reflectance for the Rededge-M images were manually performed, as indicated in the previous chapter. The workflow then proceeded with the identification of the nadir pixels ( Figure 2) and the calculation of their relative sensor and sun angles ( Figure 3). For each image, the distance in pixel units between the image optical centre (provided by the manufacturer factory calibration) and the pixel coordinates representing the surface point where the sensor was at the nadir position were calculated as Equations (1) and (2): where u is the pixel column number, v the pixel row number, F the camera focal distance in pixel unit, ϕ the adjusted sensor Y orientation, ω the adjusted X sensor orientation. This allowed to locate the nadir pixel coordinates ( Figure 2). The ANIF map for each image was then calculated as the ratio between average reflectance in a 5 degree diameter around the nadir and each other pixel. The adjusted sensor azimuth κ and the solar azimuth from R solar almanac package "suncalc" [20] were then used to calculate, for each image, the map of the relative azimuth ( Figure 3) between each image pixel, the nadir pixel, and the sun azimuth (Rel_az). Relative azimuths were calculated in a 0-180 degree reference system with 0 degrees coinciding with the direction of the solar rays and 180 degrees opposing the sun. The sensor zenith angle (PS_zen) was calculated as the angle between each pixel, the focal point, and the nadir pixel by applying the law of cosines in Equation (3): from the reference system of the internal camera parameters, where FP is the distance between the focal point and each pixel, FN the distance between the focal point and the nadir, and NP the distance between the nadir pixel and each other pixel (all distances in pixel units). Finally, image maps of the relative sensor-pixel-sun zenith (Rel_zen) was then calculated as Equation (4): where Sol_zen is the solar zenith from the solar almanac, calculated from the vertical. From the maps obtained in the image calibration dataset, the ANIF values were extracted and analysed in correlation to Rel_az, Rel_zen and PS_zen values. The calibration data from the 12 selected flights were merged to represent the entire available range of hemispherical illumination and observation angles (1000 regularly distributed ANIF values per dataset, for a total of 12,000 calibration values). Based on the observed relationships, several model equations with different combinations of input angles were tested and com-Remote Sens. 2022, 14, 537 6 of 15 pared for the best ANIF prediction. The performance of the best model was then assessed in the calibration procedure and subsequently on three validation datasets (1000 regularly distributed ANIF values per dataset, for a total of 3000 calibration values), obtained in the same way as for the calibration dataset. R 2 , Root Mean Square Error (RMSE) and Mean Relative Error % (MRE) were used as model fit evaluators. The general workflow adopted for the calibration and application of the ANIF correction is presented in Figure 4.
From the maps obtained in the image calibration dataset, the ANIF values were extracted and analysed in correlation to Rel_az, Rel_zen and PS_zen values. The calibration data from the 12 selected flights were merged to represent the entire available range of hemispherical illumination and observation angles (1000 regularly distributed ANIF values per dataset, for a total of 12,000 calibration values). Based on the observed relationships, several model equations with different combinations of input angles were tested and compared for the best ANIF prediction. The performance of the best model was then assessed in the calibration procedure and subsequently on three validation datasets (1000 regularly distributed ANIF values per dataset, for a total of 3000 calibration values), obtained in the same way as for the calibration dataset. R 2 , Root Mean Square Error (RMSE) and Mean Relative Error % (MRE) were used as model fit evaluators. The general workflow adopted for the calibration and application of the ANIF correction is presented in Figure 4. . Visualization of the relative angles extracted for BRDF modelling within the sun-pixelsensor framework: (a) relative azimuth "Rel_az"; (b) sensor zenith "PS_zen", solar zenith "Sol_zen", relative zenith "Rel_zen". tracted and analysed in correlation to Rel_az, Rel_zen and PS_zen values. The calibration data from the 12 selected flights were merged to represent the entire available range of hemispherical illumination and observation angles (1000 regularly distributed ANIF values per dataset, for a total of 12,000 calibration values). Based on the observed relationships, several model equations with different combinations of input angles were tested and compared for the best ANIF prediction. The performance of the best model was then assessed in the calibration procedure and subsequently on three validation datasets (1000 regularly distributed ANIF values per dataset, for a total of 3000 calibration values), obtained in the same way as for the calibration dataset. R 2 , Root Mean Square Error (RMSE) and Mean Relative Error % (MRE) were used as model fit evaluators. The general workflow adopted for the calibration and application of the ANIF correction is presented in Figure 4. . Visualization of the relative angles extracted for BRDF modelling within the sun-pixelsensor framework: (a) relative azimuth "Rel_az"; (b) sensor zenith "PS_zen", solar zenith "Sol_zen", relative zenith "Rel_zen".

Correction Evaluation
The optimal ANIF correction model was then applied to all available datasets, where the datasets that were not used for calibration were used for validation (distinction marked in Table 1Error! Reference source not found., last column). Images were first converted to reflectance, then Rel_zen and PS_zen angles were calculated to predict and apply

Correction Evaluation
The optimal ANIF correction model was then applied to all available datasets, where the datasets that were not used for calibration were used for validation (distinction marked in Table 1, last column). Images were first converted to reflectance, then Rel_zen and PS_zen angles were calculated to predict and apply the ANIF. The corrected images were then used in PIX4Dmapper for corrected orthomosaic generation, bypassing the camera radiometric calibration and reflectance calculation routine provided in the software. Two steps of correction evaluation were then executed: the first (i) on the images of both calibration and validation datasets, the second (ii) on a subset of orthomosaics, selected visually, which showed distinct "striping" effects (indicated with an * in Table 1 names).
During step (i), to quantify the efficacy of reflectance normalization at the image level, the change in reflectance of the same pixel location coordinates observed from different contiguous images (overlapping FOV but with different shooting time, position, and orientation) was studied before and after correction application. This was accomplished by selecting a set of 220 sample pixels from the orthomosaics of the calibration dataset (12 flights) and 80 sample pixels from the validation dataset (three flights) and tracking their corresponding pixel in all their "parent" images. Between 8 and 20 parent images were identified per pixel sample. The standard deviation (sd) of the reflectance of these parent pixels across different images was then compared between original and corrected images datasets. A substantial reduction in standard deviation after correction is then assumed as indicative for an efficient ANIF correction.
In step (ii), to quantify the "striping" effect in each orthomosaic, before and after the correction, spectral data were extracted from 8-10 transects positioned perpendicularly to the stripes. In this way, the stripes were expected to appear as spectral peaks and valleys along the transects. The positions of the drone flight lines were then overlaid on the transects, and this made it possible to study the magnitude of the striping effect in relation to the flight line positions. Two indexes were then extracted: a "spectral sinuosity" index (SSin), and the correlation between the magnitude of spectral peaks/valleys and the distance from the drone flight lines. SSin was calculated as Equation (5): where N is the number of pixels along the transect, p the pixel index along the transect, r the reflectance (0-100 scale) and l the average pixel length on the transect (m). It is the ratio between the transect distance weighted by its pixel-by-pixel spectral variability and the real planimetric distance, with domain SSin ≥ 1, where the minimum value of 1 means no spectral variability across transect distance. The average SSin of all transects of a map is then calculated.
To calculate the correlation between striping and distance from the flight lines, the magnitude of the peaks/valleys was associated to the residuals of a linear interpolation of the reflectance values along each transect. The correlation was then calculated between these residuals and the horizontal distance to the closest flight line. The underlying assumption is that, in an ideal multispectral image dataset free of anisotropic effects, reflectance changes should not show correlation with distance from the image capture positions.

Striping Issue
The effects of soil reflectance anisotropy on the products of our UAS-based multispectral survey were evident both in the images and in the final products of the photogrammetric orthomosaic reconstruction. A spectral gradient in apparent correlation with the direction of incoming sunlight was evident in all the images (e.g., in Figure 5a), in accordance with previous observations on the differences between the forward-and back-scattering reflectance of soil surfaces [21][22][23][24]. With the camera tilted forward in a fixed position to constantly aim at the nadir during cruising speed, this spectral gradient was invariant from the changing drone heading. We therefore excluded any other major systematic effect on measured reflectance (i.e., due to the camera orientation) than illumination direction. The combination of this systematic gradient behaviour and the regular geometry of the flight lines used for the image acquisition resulted in spectral artefacts in the final products of the photogrammetric orthomosaic reconstruction. These artefacts (example presented in Figure 5b) manifested as stripes with different spectral intensities, parallel and in concomitance with flight lines. When observed along transects perpendicular to the flight lines (e.g., in Figure 6a), these stripes took the form of spectral peaks and valleys (e.g., in Figure 6b).
The effects of soil reflectance anisotropy on the products of our UAS-based multispectral survey were evident both in the images and in the final products of the photogrammetric orthomosaic reconstruction. A spectral gradient in apparent correlation with the direction of incoming sunlight was evident in all the images (e.g., in Figure 5a), in accordance with previous observations on the differences between the forward-and backscattering reflectance of soil surfaces [21][22][23][24]. With the camera tilted forward in a fixed position to constantly aim at the nadir during cruising speed, this spectral gradient was invariant from the changing drone heading. We therefore excluded any other major systematic effect on measured reflectance (i.e., due to the camera orientation) than illumination direction. The combination of this systematic gradient behaviour and the regular geometry of the flight lines used for the image acquisition resulted in spectral artefacts in the final products of the photogrammetric orthomosaic reconstruction. These artefacts (example presented in Figure 5b) manifested as stripes with different spectral intensities, parallel and in concomitance with flight lines. When observed along transects perpendicular to the flight lines (e.g., in Figure 6a), these stripes took the form of spectral peaks and valleys (e.g., in Figure 6b).

Model Fit
The results of the calibration data extraction are summarized in Figure 7Error! Reference source not found., where the correlation plots show that both Rel_az and Rel_zen exhibited a stable behaviour with the ANIF across all the sampled conditions of solar il-

Model Fit
The results of the calibration data extraction are summarized in Figure 7, where the correlation plots show that both Rel_az and Rel_zen exhibited a stable behaviour with the ANIF across all the sampled conditions of solar illumination. Different equation arrangements of the Rel_az, Rel_zen, and PS_zen angles were tested for their overall ANIF prediction capability and for their image correction capability. The chosen form of the final model was Equation (6): with a, b, c, and d being the regression coefficients. The model was chosen for its high fit (R 2 = 0.81, RMSE = 0.07, MRE = 0% in Figure 8a) and after the assessment of the absence of artefacts in the corrected images, and of the "flattening" of the spectral value across the image field. The spectral "flattening" was visually examined through spectral data extraction in transects, as shown in Figure 9, in a subset of two images per calibration set (2 × 12). The model was tested on the three validation datasets, providing the results summarized in Figure 8b. The validation results showed that the model seems to provide stable and consistent ANIF predictions across the tested observation/illumination range (0.81 < R 2 < 0.91, 0.08 < RMSE < 0.1) but with a dataset-dependent bias, showed by a noticeable positive or negative MRE (%) and evident in the observed vs. predicted plot (Figure 8b).

Model Fit
The results of the calibration data extraction are summarized in Figure 7Error! Reference source not found., where the correlation plots show that both Rel_az and Rel_zen exhibited a stable behaviour with the ANIF across all the sampled conditions of solar illumination. Different equation arrangements of the Rel_az, Rel_zen, and PS_zen angles were tested for their overall ANIF prediction capability and for their image correction capability. The chosen form of the final model was Equation (6): ANIF ~ a•Rel_zen 1 + b•Rel_zen 2 + c•Rel_zen 3 + d•PS_zen (6) with a, b, c, and d being the regression coefficients. The model was chosen for its high fit (R 2 = 0.81, RMSE = 0.07, MRE = 0% in Figure 8a) and after the assessment of the absence of artefacts in the corrected images, and of the "flattening" of the spectral value across the image field. The spectral "flattening" was visually examined through spectral data extraction in transects, as shown in Figure 9Error! Reference source not found., in a subset of two images per calibration set (2 × 12). The model was tested on the three validation datasets, providing the results summarized in Figure 8Error! Reference source not found.b.
The validation results showed that the model seems to provide stable and consistent ANIF predictions across the tested observation/illumination range (0.81 < R 2 < 0.91, 0.08 < RMSE < 0.1) but with a dataset-dependent bias, showed by a noticeable positive or negative MRE (%) and evident in the observed vs. predicted plot (Figure 8b).

Correction Evaluation Results
Step (i) of our analysis showed that, on average, the measurement of the same pixel location from different images was more spectrally stable after the proposed ANIF correction. In the calibration dataset, the standard deviation of the measured sample pixel reflectance from the corrected images (mean = 0.023) was indeed lower (t [429] = 6.03, p < 0.05) than in the original images (mean = 0.032). This is also supported by a comparison of the density plots before and after correction ( Figure 10). A comparable result was also obtained with the validation dataset (density plots in Figure 10), were the sd of reflectance in the corrected images (mean = 0.030) was significantly lower (t [172] = 3.58, p < 0.05) than with the original images (mean = 0.037). The reflectance stability of the validation set was generally lower (corrected images obtained the same mean and sd density distribution as the original images in the calibration dataset) but the starting conditions were also worse (starting from an sd of 0.037).          in Figure 10), were the sd of reflectance in the corrected images (mean = 0.0 icantly lower (t [172] = 3.58, p < 0.05) than with the original images (mean reflectance stability of the validation set was generally lower (corrected im the same mean and sd density distribution as the original images in the taset) but the starting conditions were also worse (starting from an sd of 0. Regarding the evaluation step (ii), the average SSin index calculated f lected sets ( Figure 11) decreased in all cases except one ("Sicy"), showing proposed correction, the sinuosity is generally reduced by a minimum of 0 imum of 5.97%. A visual SSin reference example can be found in Figure  original and corrected transect provided SSin values of 1.030 and 1.022, resp a decrease in sinuosity of 0.77%. Without a clear reference sinuosity index t between a map affected by stripes and a corrected map, apart from a visual an analysis of the spectral data in relation to the drone flight lines was carr port the positive effects of our correction methodology. The correlation bet nitude of the spectral peaks and valleys on the transect and the distance flight line (sensor position) is provided in Figure 13 for the nine selected d ally, a higher correlation (both positive and negative) was observed in the where the peaks and/or valleys of the stripes seem to be either very close ample in Figure 12Error! Reference source not found.) or very far from t On the other hand, a weaker or zero correlation was generally observed on maps. The position of the sensor in the presence of the stripes did not sh behaviour: in some cases, both peaks and valleys were associated with hig the flight lines, resulting in an overall negative correlation; in other cases, only valleys were close to the flight lines, resulting in either a positive or n lation, depending on the local transect linear interpolation. In general, a low associated with the corrected orthomosaics indicated that the striping ef tively eliminated in most cases and at least reduced in some others. As a v Regarding the evaluation step (ii), the average SSin index calculated for the nine selected sets (Figure 11) decreased in all cases except one ("Sicy"), showing that with the proposed correction, the sinuosity is generally reduced by a minimum of 0.63% to a maximum of 5.97%. A visual SSin reference example can be found in Figure 12, where the original and corrected transect provided SSin values of 1.030 and 1.022, respectively, with a decrease in sinuosity of 0.77%. Without a clear reference sinuosity index to discriminate between a map affected by stripes and a corrected map, apart from a visual interpretation, an analysis of the spectral data in relation to the drone flight lines was carried out to support the positive effects of our correction methodology. The correlation between the magnitude of the spectral peaks and valleys on the transect and the distance to the closest flight line (sensor position) is provided in Figure 13 for the nine selected datasets. Generally, a higher correlation (both positive and negative) was observed in the original maps, where the peaks and/or valleys of the stripes seem to be either very close (as per the example in Figure 12) or very far from the flight lines. On the other hand, a weaker or zero correlation was generally observed on the corrected maps. The position of the sensor in the presence of the stripes did not show consistent behaviour: in some cases, both peaks and valleys were associated with high proximity to the flight lines, resulting in an overall negative correlation; in other cases, only peaks or only valleys were close to the flight lines, resulting in either a positive or negative correlation, depending on the local transect linear interpolation. In general, a lower correlation associated with the corrected orthomosaics indicated that the striping effect was effectively eliminated in most cases and at least reduced in some others. As a visual example of the difference between the original and the corrected spectral orthomosaic (e.g., rededge band) and of the elimination of the "striping" effect, the orthomosaic product of the "Hostellerie" dataset is reported in Figure 14.
of the difference between the original and the corrected spectral orthomosaic (e.g., rededge band) and of the elimination of the "striping" effect, the orthomosaic product of the "Hostellerie" dataset is reported in Figure 14.   of the difference between the original and the corrected spectral orthomosaic (e.g., rededge band) and of the elimination of the "striping" effect, the orthomosaic product of the "Hostellerie" dataset is reported in Figure 14.

Discussion
The results of analysis step (i) showed that the stability of the spectral measurement of the same pixel location across different images sensibly improved after the application of the ANIF correction. The relative improvement in sd ( Figure 10) was almost the same for the calibration and validation datasets, although the overall spectral stability was lower in the validation dataset. This was probably due to the non-optimal terrain conditions that prevented the initial choice of the three validation datasets for calibration purposes. Contrary to the calibration dataset, the fields used for the validation data were characterized by non-homogeneous soil colours, soil crusting (accumulation on the surface of finer particles with high reflectance), vegetation residue and varying slope/aspect. In this regard, the effects of the soil colour and vegetation residuals could probably be encompassed by a richer calibration sampling. On the other hand, surface slope and aspect play much more important roles in illumination geometry, which are generally tackled with topographic correction methods [25,26], but their effect on soil anisotropy was not accounted for in the context of this research. The specific ANIF model equation in this study was developed using only relative observation/illumination angles, which are slope/aspect-invariant and only dependent on sensor and sun orientation. This represents

Discussion
The results of analysis step (i) showed that the stability of the spectral measurement of the same pixel location across different images sensibly improved after the application of the ANIF correction. The relative improvement in sd ( Figure 10) was almost the same for the calibration and validation datasets, although the overall spectral stability was lower in the validation dataset. This was probably due to the non-optimal terrain conditions that prevented the initial choice of the three validation datasets for calibration purposes. Contrary to the calibration dataset, the fields used for the validation data were characterized by non-homogeneous soil colours, soil crusting (accumulation on the surface of finer particles with high reflectance), vegetation residue and varying slope/aspect. In this regard, the effects of the soil colour and vegetation residuals could probably be encompassed by a richer calibration sampling. On the other hand, surface slope and aspect play much more important roles in illumination geometry, which are generally tackled with topographic correction methods [25,26], but their effect on soil anisotropy was not accounted for in the context of this research. The specific ANIF model equation in this study was developed using only relative observation/illumination angles, which are slope/aspect-invariant and only dependent on sensor and sun orientation. This represents a limitation and calls for further development of the proposed model. Moreover, the surface roughness and the partial or total loss of radiometric signature caused by its shadows represent additional issues in need of further scrutiny [27][28][29].
The two correction evaluation indicators adopted in step (ii) seem to show that, in general, the specific "striping" artefact affecting the orthomosaics of the available dataset was greatly reduced or eliminated following the proposed correction methodology. Some uncertainties arose, such as the observation of the SSin analysis results (Figure 13), which revealed the anomaly of higher SSin values in the "Sicy" dataset compared to the other fields. For this specific case, it must be noted that the dataset was not used in model calibration as its images were not sufficiently homogeneous due to the prominent presence of tractor tracks in most of the images. Therefore, the interpretation of step (ii) indexes was confined to this specific case study, as no other studies were found to tackle this specific problem, and no proper comparison with other indicators could be drawn. Moreover, we could not present completely objective evaluations of the correctness of the analysed spectral orthomosaics, as no independent data with reference methods were available.

Conclusions
This study provided an empirical correction methodology for UAS-based multispectral imaging datasets of bare soils to account for the effects of the BRDF, based on the application of an anisotropy factor (ANIF) correction. The dataset of the images available for this study was sufficient to calibrate a simple model relating the ANIF to a wide range of sensor-pixel-illumination conditions by using only two angles: relative sensor-pixelsun zenith and relative sensor-pixel-sun azimuth. The proposed correction was tested on 15 available datasets. As a result, the employment of ANIF-corrected images for multispectral orthomosaic generation with photogrammetric software provided spectral maps free of anisotropic-related artefacts in most cases, as assessed by several ad hoc indexes. The validation with an independent image dataset, however, showed the presence of some systematic errors, the causes of which require further investigation. Plausible future directions for model enhancement include the implementation of a topographic index to account for the effects of slope and aspect, the consideration of soil roughness and its related shadows, and/or an enrichment of the model calibration with different soil conditions (vegetation residue, humidity, crusting).