Using Multispectral Drone Imagery for Spatially Explicit Modeling of Wave Attenuation through a Salt Marsh Meadow

: O ﬀ ering remarkable biodiversity, coastal salt marshes also provide a wide variety of ecosystem services: cultural services (leisure, tourist amenities), supply services (crop production, pastoralism) and regulation services including carbon sequestration and natural protection against coastal erosion and inundation. The consideration of this coastal protection ecosystem service takes part in a renewed vision of coastal risk management and especially marine ﬂooding, with an emerging focus on “nature-based solutions.” Through this work, using remote-sensing methods, we propose a novel drone-based spatial modeling methodology of the salt marsh hydrodynamic attenuation at very high spatial resolution (VHSR). This indirect modeling is based on in situ measurements of signiﬁcant wave heights (Hm 0 ) that constitute the ground truth, as well as spectral and topographical predictors from VHSR multispectral drone imagery. By using simple and multiple linear regressions, we identify the contribution of predictors, taken individually, and jointly. The best individual drone-based predictor is the green waveband. Dealing with the addition of individual predictors to the red-green-blue (RGB) model, the highest gain is observed with the red edge waveband, followed by the near-infrared, then the digital surface model. The best full combination is the RGB enhanced by the red edge and the normalized di ﬀ erence vegetation index (coe ﬃ cient of determination (R 2 ): 0.85, root mean square error (RMSE): 0.20% / m).

Due to the complexity of salt marsh meadows, which are rarely mono-specifics along European coastlines but more similar to a patchwork of species, the calculation of attenuation values specifically for each species remains limited for the global understanding of the attenuation dynamics. Furthermore, unlike studies in the literature that report attenuation values calculated only along one or more cross-shore transects distributed on the study site [11], a methodology based on remote sensing could bring an overview to evaluate the attenuation induced by a whole meadow. Considering the high reflectance of vegetation in the red edge (RE) and near-infrared (NIR) electromagnetic spectrum, the use of sensors provided with such wavebands and related index is well-suited to the need of an overview of the wave attenuation induced by a meadow. The use of these infrared wavebands makes it possible to discriminate plant functional traits, density [12][13][14], and species identification at very high spatial resolution (VHSR) using remote-sensing methods [8].
The use of drone technology for the imagery acquisition and for the generation of digital surface models (DSMs), using structure from motion (SfM) photogrammetry [15], is of many interests, especially for the coastal monitoring [16]. It allows for creating orthophotomosaics and DSMs at VHSR and very high temporal resolution (VHTR) [17]. This is highly relevant to monitoring the changes in plant communities depending on the seasonality or after a punctual event [8,18]. The monitoring of these changes is crucial owing to the temporal impact in the services that can be furnished by the ecosystems. The modeling method of wave attenuation presented in this paper is an indirect modeling methodology based on the spectral signatures of the vegetation which are a proxy of the vegetation volume that governs the water movement.
The aim of this paper is to quantify the contribution of the drone spectral bands, derived spectral index, and DSM for the indirect spatial modeling of wave attenuation. The working hypothesis was that adding infrared information (namely, RE, NIR, and normalized difference vegetation index, NDVI) to the basic red-green-blue (RGB) information can significantly improve the modeling performance. The paper is structured as follows: (1) data acquisition and processing (wave measurements and unmanned aerial system (UAS) multispectral imagery); (2) modeling the contribution of the salt marsh to wave attenuation; (3) quantify and discuss the contribution of each spectro-spatial band and band combination to the modeling.

Study Site
The bay of Mont-Saint-Michel (France) is enshrined in the Brittany-Normandy Gulf between the Cotentin Peninsula and the northern coast of Brittany called the "Emerald Coast". This bay, subjected to a strong tidal regime, belongs to the top six areas hosting the world's highest tide [19]. The bay is delimited in the south from the lowland of the Dol Marsh polder by the Duchess Anne's dike, at the foot of which thrive important salt marsh surfaces contributing to coastal protection as natural barriers against waves. These extended salt marsh meadows cover an area of approximately 40 km 2 in the bay [20]. This study is focusing on a well-vegetated 1 km 2 meadow of the western part of the bay, where hydrodynamic conditions are calmer than on the estuarine domain of the eastern part. This area forms an outgrowth compared to its direct vicinity, because of the sediment supplies carried by one of the draining channels of the Dol Marsh, just on its eastern side ( Figure 1). The study area is composed of four principal plant communities, from north to south, (1) pioneer vegetation of Spartina townsendii/Sueda maritima/Salicornia europea; (2) a short lawn of Puccinellia maritima with some Halimione portulacoides bush and Aster tripolium; (3) a semi-woody formation of Halimione portulacoides; and (4) a dense meadow dominated by Festuca rubra [21] (Table 1).

Short lawn
Puccinellia maritima

Spartina townsendii
Drones 2020, 4, x FOR PEER REVIEW 3 of 12        Measurements of significant wave heights were made from 22 to 23 January 2019, during four successive tides, using 12 mini-pressure sensors (NKE SP2T10) positioned along two transects at the study site ( Figure 2). They were fixed to an iron rod driven into the intertidal substrate. The instruments recorded tide-and wave-induced pressure with a burst duration of 9 min every 15 min, at a frequency of acquisition of 2 Hz. The 9-min sampling frequency was chosen as a compromise between a spectrum large enough to be representative, but sufficiently short, to assure a good degree of stationarity in the megatidal environment. Wave characteristics were obtained from the measured time series by spectral analysis using fast Fourier transforms. The Fourier coefficients of the free surface elevation fluctuations were obtained from corresponding coefficients computed from the pressure time series using the frequency-dependent transfer function inferred from linear theory.
These Hm0 values are enabled to calculate 40 attenuation values (N = 40), which embodied the response within the various statistical models further tested for hydrodynamic attenuation modeling.

Significant Wave Height Measurement (Hm0) as the Model Response
Measurements of significant wave heights were made from 22 to 23 January 2019, during four successive tides, using 12 mini-pressure sensors (NKE SP2T10) positioned along two transects at the study site ( Figure 2). They were fixed to an iron rod driven into the intertidal substrate. The instruments recorded tide-and wave-induced pressure with a burst duration of 9 min every 15 min, at a frequency of acquisition of 2 Hz. The 9-min sampling frequency was chosen as a compromise between a spectrum large enough to be representative, but sufficiently short, to assure a good degree of stationarity in the megatidal environment. Wave characteristics were obtained from the measured time series by spectral analysis using fast Fourier transforms. The Fourier coefficients of the free surface elevation fluctuations were obtained from corresponding coefficients computed from the pressure time series using the frequency-dependent transfer function inferred from linear theory.
These Hm0 values are enabled to calculate 40 attenuation values (N = 40), which embodied the response within the various statistical models further tested for hydrodynamic attenuation modeling.

Semi-woody formation
Halimione portulacoides Measurements of significant wave heights were made from 22 to 23 January 2019, during four successive tides, using 12 mini-pressure sensors (NKE SP2T10) positioned along two transects at the study site ( Figure 2). They were fixed to an iron rod driven into the intertidal substrate. The instruments recorded tide-and wave-induced pressure with a burst duration of 9 min every 15 min, at a frequency of acquisition of 2 Hz. The 9-min sampling frequency was chosen as a compromise between a spectrum large enough to be representative, but sufficiently short, to assure a good degree of stationarity in the megatidal environment. Wave characteristics were obtained from the measured time series by spectral analysis using fast Fourier transforms. The Fourier coefficients of the free surface elevation fluctuations were obtained from corresponding coefficients computed from the pressure time series using the frequency-dependent transfer function inferred from linear theory.
These Hm0 values are enabled to calculate 40 attenuation values (N = 40), which embodied the response within the various statistical models further tested for hydrodynamic attenuation modeling.

Significant Wave Height Measurement (Hm0) as the Model Response
Measurements of significant wave heights were made from 22 to 23 January 2019, during four successive tides, using 12 mini-pressure sensors (NKE SP2T10) positioned along two transects at the study site ( Figure 2). They were fixed to an iron rod driven into the intertidal substrate. The instruments recorded tide-and wave-induced pressure with a burst duration of 9 min every 15 min, at a frequency of acquisition of 2 Hz. The 9-min sampling frequency was chosen as a compromise between a spectrum large enough to be representative, but sufficiently short, to assure a good degree of stationarity in the megatidal environment. Wave characteristics were obtained from the measured time series by spectral analysis using fast Fourier transforms. The Fourier coefficients of the free surface elevation fluctuations were obtained from corresponding coefficients computed from the pressure time series using the frequency-dependent transfer function inferred from linear theory.
These Hm0 values are enabled to calculate 40 attenuation values (N = 40), which embodied the response within the various statistical models further tested for hydrodynamic attenuation modeling.

Significant Wave Height Measurement (Hm 0 ) as the Model Response
Measurements of significant wave heights were made from 22 to 23 January 2019, during four successive tides, using 12 mini-pressure sensors (NKE SP2T10) positioned along two transects at the study site ( Figure 2). They were fixed to an iron rod driven into the intertidal substrate. The instruments recorded tide-and wave-induced pressure with a burst duration of 9 min every 15 min, at a frequency of acquisition of 2 Hz. The 9-min sampling frequency was chosen as a compromise between a spectrum large enough to be representative, but sufficiently short, to assure a good degree of stationarity in the megatidal environment. Wave characteristics were obtained from the measured time series by spectral analysis using fast Fourier transforms. The Fourier coefficients of the free surface elevation fluctuations were obtained from corresponding coefficients computed from the pressure time series using the frequency-dependent transfer function inferred from linear theory.
These Hm 0 values are enabled to calculate 40 attenuation values (N = 40), which embodied the response within the various statistical models further tested for hydrodynamic attenuation modeling.

Multispectral Drone Imagery Acquisition and Processing
A drone imagery campaign was undertaken on 22 March 2019, using a Sensefly eBee+® fixedwing equipped with two optical sensors: a classical RGB sensor and a Parrot Sequoia®, exhibiting an RGB (R: 660 nm; G: 520 nm; B: 450 nm, Figure 3a, b, and c), RE (730-740 nm, Figure 3d) and NIR (770-810 nm, Figure 3e) sensors ( Table 2). The imagery collection followed a flight plan whose parameters are provided in Table 3. Leveraging global navigation satellite system (GNSS) coordinates, the images were further georeferenced to 9 ground control points (XYZ) acquired through Differential-GNSS (D-GNSS) equipment (Topcon HiPer V) with horizontal and vertical accuracy of 0.02 m. An orthophotomosaic and a DSM were then derived using the SfM photogrammetry procedure with the Pix4Dmapper® software. The accuracy reached 0.05 m, 0.04 m and 0.12 m for X, Y and Z coordinates, respectively.    Table 2). The imagery collection followed a flight plan whose parameters are provided in Table 3. Leveraging global navigation satellite system (GNSS) coordinates, the images were further georeferenced to 9 ground control points (XYZ) acquired through Differential-GNSS (D-GNSS) equipment (Topcon HiPer V) with horizontal and vertical accuracy of 0.02 m. An orthophotomosaic and a DSM were then derived using the SfM photogrammetry procedure with the Pix4Dmapper ® software. The accuracy reached 0.05 m, 0.04 m and 0.12 m for X, Y and Z coordinates, respectively. Radiometric corrections were also applied accounting for the optical instrument factors (vignetting, spectral response, etc.), the variation in solar irradiance and angle, using a radiometric target to provide top of canopy reflectance [22].

Multispectral Drone Imagery Acquisition and Processing
Drones 2020, 4, x FOR PEER REVIEW 6 of 12 Radiometric corrections were also applied accounting for the optical instrument factors (vignetting, spectral response, etc.), the variation in solar irradiance and angle, using a radiometric target to provide top of canopy reflectance [22].

Modeling the Contribution of the Salt Marsh to Wave Attenuation
Data extracted from the drone imagery were considered as descriptors of the salt marsh composition and 3D configuration, thus predictors for the statistical modeling. Using the 1D Hm0 in situ values and the 2D seven predictors (R, G, B, RE, NIR, NDVI, DSM), an array of linear regression models were achieved to spatially predict the hydrodynamic attenuation (Y) induced by the salt marsh meadow. The contribution of each predictor (βi) was therefore independently tested through single regressions (Table 4): Different combinations (Table 3) of predictors were then tested using multiple regressions:

Modeling the Contribution of the Salt Marsh to Wave Attenuation
Data extracted from the drone imagery were considered as descriptors of the salt marsh composition and 3D configuration, thus predictors for the statistical modeling. Using the 1D Hm 0 in situ values and the 2D seven predictors (R, G, B, RE, NIR, NDVI, DSM), an array of linear regression models were achieved to spatially predict the hydrodynamic attenuation (Y) induced by the salt marsh meadow. The contribution of each predictor (β i ) was therefore independently tested through single regressions (Table 4): Different combinations (Table 3) of predictors were then tested using multiple regressions: The coefficient of determination R 2 and the root mean square error value (RMSE) were used to assess the relevance of each model tested.

Results
The 23 linear regression models, both 7 simple and 16 multiple, allowed to quantify the predictors' contribution to the hydrodynamic attenuation modeling ( Table 5).
Dealing with multiple linear regressions, each predictor was added to the classic RGB combination to quantify how much each one could improve the coefficient of determination of the RGB model (R 2 : 0.54). In so doing, the combinations with the three highest R 2 were RGB + RE (R 2 : 0.73), RGB + NIR (R 2 : 0.71), RGB + DSM (R 2 : 0.64).
For this reason, the spatially explicit model was achieved using the RGB + RE + NDVI combination ( Figure 4).
The spatially explicit-model drew attenuation on values displaying important variations, from an increase of Hm 0 of 2.5%/m to a decrease of 3.5%/m on the study site. The mean Hm 0 attenuation value on the study site approximated 0.25%/m. The simple linear regressions results showed the relevance of each predictor to explain the wave attenuation, via its coefficient of determination, ranked in ascendant order as follows: RE (R 2 : 0.24), DSM (R 2 : 0.29), NIR (R 2 : 0.32), R (R 2 : 0.33), NDVI (R 2 : 0.41), B (R 2 : 0.50) and G (R 2 : 0.51).
Dealing with multiple linear regressions, each predictor was added to the classic RGB combination to quantify how much each one could improve the coefficient of determination of the RGB model (R 2 : 0.54). In so doing, the combinations with the three highest R 2 were RGB + RE (R 2 : 0.73), RGB + NIR (R 2 : 0.71), RGB + DSM (R 2 : 0.64).
For this reason, the spatially explicit model was achieved using the RGB + RE + NDVI combination (Figure 4). The spatially explicit-model drew attenuation on values displaying important variations, from an increase of Hm0 of 2.5%/m to a decrease of 3.5%/m on the study site. The mean Hm0 attenuation value on the study site approximated 0.25%/m.

Discussion
This experiment highlighted interest in the use of infrared wavebands (RE, NIR and derived data NDVI) and topographic (DSM) data in addition to RGB channels for the spatial modeling of the wave attenuation through a salt marsh meadow.
Results showed that the individual most effective predictors (Table 4), to explain the attenuation values calculated from Hm 0 , in descending order, correspond to the reflectance (G and NIR) and absorbance (R and B) peaks of chlorophyll pigments [23], their NDVI combination (vegetation density), the surface changes (DSM) and the RE.
The addition of individual predictors to the RGB combination allowed us to observe some interesting gains concerning the coefficient of determination ( Figure 5), especially the addition of RE and NIR predictors that improved by +0.19 and +0.17, respectively, the R 2 of the multiple linear regression of the RGB model. Thus, from the predictor with the highest gain to the predictor with the lowest gain, there was: RE (+0.19), NIR (+0.17), DSM (+0.10) and NDVI (+0.04). and NIR predictors that improved by +0.19 and +0.17, respectively, the R of the multiple linear regression of the RGB model. Thus, from the predictor with the highest gain to the predictor with the lowest gain, there was: RE (+0.19), NIR (+0.17), DSM (+0.10) and NDVI (+0.04).
These results can be explained by the fact that these predictors highlighted the natural vegetation elements that increase the roughness or represent topographical changes able to mitigate waves. The best performance of the RE and NIR predictors, witnessing the vegetation, indicate that the structural complexity of vegetation seems more effective to model the attenuation than the topographical change, highlighted by the DSM. The contribution of the NDVI appeared lower, perhaps due to the R integration into it, generating an information redundancy with the RGB combination.
With regard to the models based on the RGB dataset combined with multiple predictors, the best gains were observed with the addition of the RE + NDVI (+0.31), RE + DSM (+0.30), RE + NDVI + DSM (+0.29), RE + NIR + NDVI (+0.27), RE + NIR + DSM (+0.26) predictors ( Figure 5). These results highlighted the ubiquity of the RE predictor into the most efficient combination predictors for the spatially-explicit modeling of the Hm0 attenuation, because of its high sensitivity to a medium to high level of chlorophyll content, complementary to the RGB spectral signatures of the salt marsh vegetation [24].
Results of the combinations RGB + RE + NDVI (R 2 : 0.85), RGB + RE + DSM (R 2 : 0.84) and RGB + RE + NDVI + DSM (R 2 : 0.83) ( Figure 5) were strongly comparable. It implies that the NDVI and the DSM bring the same quantity of complementary information to the modeling. Once the RE intermediate reflectance is integrated, the vegetation density (NDVI) and the topographic features are almost identical for the attenuation prediction.  These results can be explained by the fact that these predictors highlighted the natural vegetation elements that increase the roughness or represent topographical changes able to mitigate waves. The best performance of the RE and NIR predictors, witnessing the vegetation, indicate that the structural complexity of vegetation seems more effective to model the attenuation than the topographical change, highlighted by the DSM. The contribution of the NDVI appeared lower, perhaps due to the R integration into it, generating an information redundancy with the RGB combination.
With regard to the models based on the RGB dataset combined with multiple predictors, the best gains were observed with the addition of the RE + NDVI (+0.31), RE + DSM (+0.30), RE + NDVI + DSM (+0.29), RE + NIR + NDVI (+0.27), RE + NIR + DSM (+0.26) predictors ( Figure 5). These results highlighted the ubiquity of the RE predictor into the most efficient combination predictors for the spatially-explicit modeling of the Hm0 attenuation, because of its high sensitivity to a medium to high level of chlorophyll content, complementary to the RGB spectral signatures of the salt marsh vegetation [24].
Results of the combinations RGB + RE + NDVI (R 2 : 0.85), RGB + RE + DSM (R 2 : 0.84) and RGB + RE + NDVI + DSM (R 2 : 0.83) ( Figure 5) were strongly comparable. It implies that the NDVI and the DSM bring the same quantity of complementary information to the modeling. Once the RE intermediate reflectance is integrated, the vegetation density (NDVI) and the topographic features are almost identical for the attenuation prediction.
The highest attenuation values (between 1.5 and 3.5%/m) were mainly located in two areas of the study site: (1) on the marsh edge, probably due the abrupt topographical change of the slope, which induced a wave overwash, and (2) on the high marsh where the vegetation was the tallest on the study site. By contrast, the lowest attenuation values (between −2.5 and 0.5%/m) were observed where the vegetation density was lower or quasi absent (i.e., mudflat), and in the channels in front of the marsh where the bathymetry was higher ( Figure 6).
The highest attenuation values (between 1.5 and 3.5%/m) were mainly located in two areas of the study site: 1) on the marsh edge, probably due the abrupt topographical change of the slope, which induced a wave overwash, and 2) on the high marsh where the vegetation was the tallest on the study site. By contrast, the lowest attenuation values (between -2.5 and 0.5%/m) were observed where the vegetation density was lower or quasi absent (i.e., mudflat), and in the channels in front of the marsh where the bathymetry was higher ( Figure 6).
The results indicate also that vegetation is (through its standing biomass) a better proxy than DSM. It can be explained by the dense standing biomass (with a specific spectral signature), which translates into greater roughness than the DSM information, probably underestimating the structural complexity of the biomass because of the dense canopy closure of this low vegetation.
The mean Hm0 attenuation value on the study site is comparable to the values obtained through north-west European mixed salt marshes (0.34%/m [3], 0.25-0.30%/m [25]). The slightly lower results can be explained by the study area that presents a long mudflat, reducing the mean attenuation value on the whole site.
Some artefacts due to a vignetting effect following the flight lines appeared on the spatially explicit model. This noise can be nevertheless retrieved from the original B, G, R, RE and NIR data sources using a principal component analysis (PCA), for further investigation (Figures 6). This vignetting effect could be imputed to the sensors and especially the Sequoia® sensor which exhibits some troubles into the radiometric correction, and to the lateral overlap ratio. Concerning the future improvements of this modeling methodology, the Hm0 data acquisition could be enhanced by increasing the number of measurement stations and by widening their distribution on the study site. Augmentation of the side overlap ratio, from 65% to 80%, and the use of the second-generation sensor (e.g., Parrot Sequoia+®, leveraging an automatic radiometric calibration system with a shorter shutter interval) is advised, to reduce the vignetting effect. These corrections hold great potential to ameliorate the prediction of the built regression models.
Complementary topographical predictors could also be added to the models such as the roughness index or slope index sourced from the photogrammetry methods, or from other sensors like airborne light detection and ranging system (LiDAR), bringing information about the vegetation topographical but also spectral (intensity) signatures [26,27]. The spatially explicit modeling could also be enhanced using artificial or convolutional neural networks instead of linear regressions.

Conclusions
The prediction of multispectral drone imagery was very efficient and accurate in achieving spatially explicit modeling of the wave attenuation in a salt marsh meadow, in addition to in situ wave measurements. Individual G, B, NDVI, R, NIR, DSM and finally RE bands decreasingly brought insights into the attenuation simple linear modeling, emphasizing the relevance of the reflectance The results indicate also that vegetation is (through its standing biomass) a better proxy than DSM. It can be explained by the dense standing biomass (with a specific spectral signature), which translates into greater roughness than the DSM information, probably underestimating the structural complexity of the biomass because of the dense canopy closure of this low vegetation.
The mean Hm 0 attenuation value on the study site is comparable to the values obtained through north-west European mixed salt marshes (0.34%/m [3], 0.25-0.30%/m [25]). The slightly lower results can be explained by the study area that presents a long mudflat, reducing the mean attenuation value on the whole site.
Some artefacts due to a vignetting effect following the flight lines appeared on the spatially explicit model. This noise can be nevertheless retrieved from the original B, G, R, RE and NIR data sources using a principal component analysis (PCA), for further investigation ( Figure 6). This vignetting effect could be imputed to the sensors and especially the Sequoia ® sensor which exhibits some troubles into the radiometric correction, and to the lateral overlap ratio.
Concerning the future improvements of this modeling methodology, the Hm 0 data acquisition could be enhanced by increasing the number of measurement stations and by widening their distribution on the study site. Augmentation of the side overlap ratio, from 65% to 80%, and the use of the second-generation sensor (e.g., Parrot Sequoia+ ® , leveraging an automatic radiometric calibration system with a shorter shutter interval) is advised, to reduce the vignetting effect. These corrections hold great potential to ameliorate the prediction of the built regression models.
Complementary topographical predictors could also be added to the models such as the roughness index or slope index sourced from the photogrammetry methods, or from other sensors like airborne light detection and ranging system (LiDAR), bringing information about the vegetation topographical but also spectral (intensity) signatures [26,27]. The spatially explicit modeling could also be enhanced using artificial or convolutional neural networks instead of linear regressions.

Conclusions
The prediction of multispectral drone imagery was very efficient and accurate in achieving spatially explicit modeling of the wave attenuation in a salt marsh meadow, in addition to in situ wave measurements. Individual G, B, NDVI, R, NIR, DSM and finally RE bands decreasingly brought insights into the attenuation simple linear modeling, emphasizing the relevance of the reflectance and absorbance peaks, the density of the vegetation communities, and then the induced surface changes. Combined with the RGB standard dataset (R 2 : 0.54), the highest gain by individual predictors was reached by the RE (+0. 19), followed by the NIR (+0.17), then the DSM (+0.10). The joint full-combinations highlighted the presence of the RE into the 5 best datasets, with the highest one composed of the RGB + RE + NDVI (R 2 : 0.85), closely followed by the RGB + RE + DSM (R 2 : 0.84) and RGB + RE + NDVI + DSM (R 2 : 0.83). These last results indicated a relative redundancy in the information drawn from the NDVI and the DSM.
In this experimental model, significant wave height values ranged from an increase of 2.5%/m to a decline of 3.5%/m. The modeling allowed us to identify what and where are the most effective areas to attenuate the waves, providing strong advice for future wetland restoration programs.