Next Article in Journal
Rapid Flood Mapping and Evaluation with a Supervised Classifier and Change Detection in Shouguang Using Sentinel-1 SAR and Sentinel-2 Optical Data
Previous Article in Journal
Using Linear Regression, Random Forests, and Support Vector Machine with Unmanned Aerial Vehicle Multispectral Images to Predict Canopy Nitrogen Weight in Corn
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mapping Benthic Habitats by Extending Non-Negative Matrix Factorization to Address the Water Column and Seabed Adjacency Effects

1
CNRS, Centrale Marseille, Institut Fresnel, Aix Marseille Univ, F-13013 Marseille, France
2
SeaTech, LIS Laboratory, Université de Toulon, CNRS UMR 7020, 83041 Toulon, France
3
IRAP, Observatoire Midi-Pyrénées, Université de Toulouse, UPS-CNRS-OMP, 14 Av. Edouard Belin, 31400 Toulouse, France
4
Laboratoire Atmosphères Milieux Observations Spatiales (LATMOS), Sorbonne Université, CNRS-INSU, Boulevard de L’Observatoire, CS 34229, CEDEX 06304 Nice, France
5
CS Systemes d’Information Company, 31506 Toulouse CEDEX 05, France
6
ONERA/DOTA, BP 74025-2 Avenue Edouard Belin, FR-31055 Toulouse CEDEX 4, France
7
DGA/DS/MRIS, 75509 Paris CEDEX 15, France
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(13), 2072; https://doi.org/10.3390/rs12132072
Submission received: 4 June 2020 / Revised: 19 June 2020 / Accepted: 23 June 2020 / Published: 27 June 2020
(This article belongs to the Section Remote Sensing Image Processing)

Abstract

:
Monitoring of coastal areas by remote sensing is an important issue. The interest of using an unmixing method to determine the seabed composition from hyperspectral aerial images of coastal areas is investigated. Unmixing provides both seabed abundances and endmember reflectances. A sub-surface mixing model is presented, based on a recently proposed oceanic radiative transfer model that accounts for seabed adjacency effects in the water column. Two original non-negative matrix factorization ( N M F )-based unmixing algorithms, referred to as W A D J U M (Water ADJacency UnMixing) and W U M (Water UnMixing, no adjacency effects) are developed, assuming as known the water column bio-optical properties. Simulations show that W A D J U M algorithm achieves performance close to that of the N M F -based unmixing of the seabed without any water column, up to 10 m depth. W U M performance is lower and decreases with the depth. The robustness of the algorithms when using erroneous information about the water column bio-optical properties is evaluated. The results show that the abundance estimation is more reliable using W A D J U M approach. W A D J U M is applied to real data acquired along the French coast; the derived abundance maps of the benthic habitats are discussed and compared to the maps obtained using a fixed spectral library and a least-square ( L S ) estimation of the seabed mixing coefficients. The results show the relevance of the W A D J U M algorithm for the local analysis of the benthic habitats.

Graphical Abstract

1. Introduction

Remote measurement of benthic reflectance addresses a range of compelling environmental challenges, such as monitoring coral reefs [1], seagrass [2], or invasive seaweed [3], which serve key roles for coastal populations and the global environment. Coastal and inland waters are generally complex environments, whose remotely sensed reflectance may be highly variable due to simultaneous changes in bathymetry, water quality, seabed type, water surface and atmospheric conditions. In shallow waters, the decoupling of these effects has been shown to be more accurate when using hyperspectral data instead of multispectral data [4,5,6]. This is because a higher number of spectral bands, as well as an increased spectral resolution, allow reducing confounding effects between optically active parameters, e.g., by detecting the subtle changes in reflectance that originate from narrow absorption regions potentially present in seabed albedo [7,8]. Therefore, hyperspectral airborne remote sensing has been widely used for mapping water composition and bathymetry [9,10]. In coastal environments, hyperspectral remote-sensing methods that allow the simultaneous retrieval of bathymetry, water quality and benthic cover are often based on the inversion of a radiative transfer model that describes how light propagates into the water [11]. Semi-analytic parametric models [12,13], as well as numerical models such as Hydrolight [14], or more recently, the O S O A A open source radiative transfer model [15], have been developed to describe the relationship between the seabed reflectance and remote-sensing reflectance by taking into account water attenuation. These models usually consider that four parameters affect the water-leaving radiance in addition to the water molecules: seabed depth and concentrations of three optically active constituents, namely chlorophyll ( C h l ), suspended non-algal particles ( S M ), and colored dissolved organic matter ( C D O M ) [16]. A classical approach could rely on estimating the parameters of the model, which provide the water column constituents (concentrations of chlorophyll, suspended matter, dissolved organic matter , and bathymetry). The retrieval of the bio-optical properties of the water column using an inversion algorithm is generally obtained using either look-up tables, iterative least-square (LS) [16], Maximum Likelihood (ML) [17,18], or Maximum A Posteriori (MAP) optimization [19,20].
In shallow water, the water column constituents estimation requires knowledge of the seabed reflectance spectrum. Generally, due to the spatial resolution of airborne or space sensors and to the complexity and variability of coastal areas, the seabed pixel reflectance is supposed to be a linear combination of endmember spectra.
When the seabed reflectance is not known but the endmember spectra are available, the coefficients of the linear mixing can be additionally estimated through the inversion process [21]. In the case when a sum-to-unity constraint is applied to the coefficients of the linear mixture, they represent the percentages of each spectrum in the considered pixel reflectance, also called abundance fractions (or abundances). This makes it possible to obtain abundance maps which provide understandable and clear information about the seabed composition. When only the shape of the endmember spectra is available, unknown scale factors are also necessary to provide the true endmember spectra. Then the mixing coefficients may depend both on those unknown scale factors as well as on the proportions of the materials, so the sum-to-unity constraint should not be applied, the estimation thus providing only coefficients maps for which the interpretation is not straightforward. Although the spectral library might not represent properly all the pure materials, it could be sufficiently comprehensive to be used as a basis. Then, the mixing coefficients may still be used to estimate the seabed spectrum. Note, however, that the coefficients derived for the latter case do not provide maps that are representative of the benthic species distribution.
When a spectral library is not available or not sufficiently representative of the benthic species, it is necessary to estimate both the mixing coefficients, and additionally the endmember signatures. Here, a difference is made between unsupervised unmixing, which we refer to estimating both mixing coefficients (abundances if the sum-to-unity constraint is applied), and endmember spectra in each observed pixel, and coefficient estimation, which we refer to estimating only the mixing coefficients, using a spectral library. Actually, due to natural variability, existing spectral libraries may not realistically represent the pure materials for a given observation, thus leading to erroneous abundance estimation. Interestingly, the non-negative matrix factorization ( N M F ) class of methods, broadcasted by Lee and Seung in the seminal articles [22,23], and extensively developed in [24], makes it possible to estimate both abundances and endmembers. Furthermore, such a method has shown to be able to adapt to various mixing models [25,26]. As an example, the case of multiple reflections between endmembers has been recently addressed using N M F [27,28,29]. The consideration of the inter-pixel variability also remains a challenge for the community, and it has been addressed using unmixing methods [30,31,32].
In the case of seabed analysis, the influence of the water column needs to be included in the mixing model. Generally, the water constituents and the bathymetry are estimated. In [33], J. Goodman et al. introduced mixing coefficient estimation for seabed, using an adapted version of Lee’s parametric semi-analytical model [12,21,34], and linear mixing for seabed reflectance, using a known spectral library. They successively perform Lee’s inversion and estimate the seabed spectra mixing coefficients, in an algorithm called L I G U . A variant of L I G U is also used in the Z y C I A T toolbox by [35]. Seabed mapping is also developed in [36]. In [37], Maritorena’s radiative transfer model [13] has been used to combine iteratively M L estimation and N M F unmixing to obtain simultaneously the water parameters, the endmember spectra and the relative abundances. Recently a Bayesian approach has been proposed [20] for linear mixtures to jointly estimate the seabed reflectance and water optical properties while flexibly incorporating varied domain knowledge and in situ measurements, using Maritorena’s model, and simulations from Hydrolight radiative transfer code. A brief review of seabed analysis methods can be found in [38].
The adjacency effects represent the contamination of the radiance of a given target by photons scattered by the surrounding environment of this target. Most of the previous studies dealing with the adjacency effects in ocean color remote sensing, aimed at analyzing the influences of the atmosphere and/or of the neighboring land pixels on the radiance measured by a satellite sensor [39,40,41,42]. In [43], the authors propose an Alternating Direction Method of Multipliers ( A D M M )-based unmixing method that takes into account adjacency effects in the atmosphere. Authors [44] have investigated the detection of adjacency effects in the ground implying water areas, as being a general unmixing problem, using several N M F -based methods. In [45], the focus is made for the first time on the analysis of the influence of the seabed adjacency effects; an analytical formulation of the radiative transfer equation is developed that includes the adjacency effects. The radiative transfer model developed in [45] has been used in a previous work [46], to introduce the new specific N M F -based unmixing method that is developed here.
The objectives of this paper are (i) to propose two sub-surface mixing models, either accounting for the adjacency effects or not, (ii) to develop two adapted N M F -based unmixing methods, in accordance with these models, iii) to analyze and quantify the efficiency of the proposed methods, by evaluating the impact of the adjacency effects on the results, and iv) to obtain abundance maps of the seabed from airborne hyperspectral images of a coastal area and assess their relevance.
To reach these objectives, the analytical formulation of the radiative transfer into the water column, recently proposed in [45], is taken into account as the original complete model, which is then adapted to the unmixing purpose. The adapted radiative transfer model and the mixed pixel decomposition are accounted for to derive a specific sub-surface mixing model. A new N M F -based unmixing method, hereafter referred to as W A D J U M (Water ADJacency UnMixing), is proposed. A simpler variant, referred to as W U M (Water UnMixing) algorithm, which does not account for the adjacency effects induced by the seabed and the hydrosols present in the water column, is also presented. An analysis of the performance of the proposed unmixing methods is carried out by assessing the relative importance of the adjacency effects, and also the relative importance of the presence of a water column. Additionally, the impact of an error made on the water column composition is studied. Last, the W A D J U M algorithm is applied to real original data acquired on the coast of the Porquerolles Island, South of France. The resulting abundance maps are evaluated by comparison with a ground truth, and with coefficient maps that are obtained using an L S coefficient estimation approach and a spectral library.
The paper is organized as follows: Section 2 is devoted to the theoretical developments. In Section 3, the unmixing performance is assessed by means of simulations that are performed using the O S O A A open source radiative transfer model [15]. Section 4 provides results of the application of the proposed unmixing methods to real data, and Section 5 is devoted to a discussion.

2. Models and Methods

In this section, the underwater radiative transfer equation that was recently developed in [45] is briefly described. Such an analytic formulation has been established from the complete equation of the radiative transfer into the water column. This previous work has highlighted the influence of the adjacency effects induced by the seabed and by the hydrosols present in the water column.
The notations used in the paper are given in Table 1.

2.1. Radiative Transfer Equation that Makes the Adjacency Effects Explicit

This work follows an in-depth analysis of various sources of photons that contribute to the sub-surface upward radiance radiative transfer phenomena within the water column, presented in [45]. From this paper, a schematic representation of the upward radiances corresponding to direct, diffuse, and water-leaving radiance is shown in Figure 1a. Let us note that the direct sun reflection on the water surface does not appear in this schematic representation, since aerial measurements are commonly acquired out of the sunglint geometry; at worse, the sunglint radiance contribution could be corrected from the data using dedicated sunglint correction algorithms [47].
The environment of the target pixel P of the seabed (for example, adjacent pixel M) contributes to the observed sub-surface radiance, due to the scattering of the light within the water column. Then, the total contributing seabed reflectance is not only associated with the target pixel, but also with its neighbors.
The observed sub-surface (depth z = 0 ) total upward radiance can be expressed as the sum of three main terms (Figure 1a):
L u ( P , 0 ) = L d i r ( P , 0 ) + L d i f ( P , 0 ) + L u d p ( P , 0 ) .
L u d p is the radiance due to the upward photons that have not interacted with the seabed. Following the previous work in [45], the upward radiances L d i r and L d i f are both related to the downward flux, to the seabed reflectance and to the water column transmittance:
L u ( P , 0 ) = E t o t ( P ) π T d i r ( 0 ) ρ t + E t o t ( P ) π T d i f ( 0 ) ρ e n v + L u d p ( P , 0 ) .
The symbol ⊙ denotes the elementwise product and each bold term corresponds to a spectral vector. ρ t and ρ e n v respectively refer to the target (pixel P) reflectance and the environment (P and adjacent pixels) reflectance.

2.2. Downward Flux E t o t

The total incident downward flux that reaches the pixel P, E t o t ( P ) , is the sum of four components, corresponding respectively to the direct, diffuse, coupling and reflected light (see Figure 1a):
E t o t = E d i r + E d i f + E c o u p + E r e f .
The coupling term E c o u p depends on multiple reflections between the seabed and the water column. It therefore depends on the reflectance of the seabed environment, including the reflectance of the target pixel, through multiple products of seabed reflectances. To evaluate the importance of this non-linear effect in reflectance, a test case has been simulated, composed of 2400 various linear mixtures of four endmember spectra that are contained in the spectral library used for the experiments (shown in Section 3, Figure 5f, dotted spectra). The mixtures are representative of the spectral variability of seabed areas in the Porquerolles Alicaster Bay (Var, France). The objective is to evaluate the amount of variability of E t o t due to the multiple reflections between the seabed and the water column, described by the coupling term E c o u p . A variability index, defined as the maximum ratio over the wavelengths between the standard deviation and the mean value of the 2400 pixels reflectance, is found as m a x r ( R ) = 61 % for the seabed. The value of E t o t has been generated using O S O A A model, including all the multiple reflections, for each pixel position and for three cases. The water quality is defined by the concentrations of the main optical constituents: chlorophyll-a, suspended matter, and colored dissolved organic matter. The clear water condition is defined as [ C h l = 0.03 mg · m 3 , S M = 0.01 mg · L 1 , C D O M = 0.01 m 1 ], and the moderately turbid is defined as [ C h l = 1 mg · m 3 , S M = 1 mg · L 1 , C D O M = 0.1 m 1 ]. The three considered cases are c 1 : clear water, seabed depth z = 1 m , c 2 : moderately turbid water, z = 1 m and c 3 : moderately turbid water, z = 10 m . For each case, the standard deviation and the mean value of E t o t over the 2400 pixels have been calculated, and the resulting variability index, m a x r ( E t o t ) , is reported in Table 2.
The coupling term has the highest influence in the case of clear water and low depth. Out of all cases considered, a variability index of 61 % for the spectra of the seabed corresponds to a maximum variability index of 1.3 % . Consequently, in the remainder of the paper, it is considered that the variation of E t o t with the seabed composition can be assimilated to (weak) noise. Then multiple reflections that cause the variations of E t o t can be neglected in the expression of the radiative transfer model of Equation (2). In the following of the paper E t o t is assumed to depend only on the water column bio-optical properties (bathymetry and turbidity), not on the seabed reflectance, so multiple products of seabed reflectances do not appear in the model.

2.3. Environment Function and Environment Reflectance Analysis

In [45], the authors make use of the environment function G e n v , as given in [42], also called δ . The environment function determines at each wavelength the weight of the target contributing to the upward diffuse radiance, relatively to its neighborhood. Typically, a value of one means that the target P contributes at 100 % to the sub-surface upward radiance; in other words, there is no adjacency effect from the neighborhood of the target for such a case. The environment function G e n v depends on the size of the target, on the optical thickness of the water layer, and on the wavelength. It is lower for smaller targets, for a high number of scattering events of hydrosols, and for deep waters. For more details the reader is encouraged to refer to [45]. Here the notation δ denotes its vectorial representation over the wavelengths, hereafter referred to as the ’environment vector’, which is plotted in Figure 1b for the case of moderately turbid water, a target pixel radius of 0.5 m, and various seabed depths. The variations of δ with the wavelength are very smooth. In consequence it is assumed in the following that the spectral vector δ can be replaced by its mean value over the wavelength domain, δ = 1 L λ = 1 L δ ( λ ) , so-called the environment parameter. Please note that δ , which provides the global contribution of the target pixel to the diffuse light in the range of 400 nm–700 nm, is approximately 72% for 5 m depth, whereas it is only 55% at 10 m, thus meaning that the adjacent pixels contribute 28% and 45% respectively to the diffuse light in these cases. Consequently, the contribution of pixels adjacent to the target pixel cannot be neglected in the diffuse sub-surface upward radiance.
Following [45], the equivalent environment reflectance of the seabed target pixel is developed as:
ρ e n v = δ ρ t + ( 1 δ ) ρ a d j ,
where ρ t is the target reflectance and ρ a d j is the reflectance of the pixels adjacent to the target. To build the underwater mixing model, it is necessary to decompose ρ t and ρ e n v as mixtures of endmember spectra.

2.4. Decomposition of Target and Environment Reflectances

In this study, the assumption of linear mixing is made for the seabed sub-pixel decomposition. This is justified by the following considerations.
(i) The linear mixing model is currently used for seabed estimation with other methods (such as L S inversion with a spectral library) [9,18,20,21,33,35,36,37,38,46].
(ii) As mentioned previously in Section 2.2, O S O A A simulations have shown that non-linear interactions such as multiple reflections between the seabed and the water column do not significantly modify the downward flux E t o t , and consequently the upward fluxes at the sea surface level. This is due to the water attenuation. At the sub-surface level, simulations showed that the only notable bottom interaction with the water column is the adjacency effect resulting from the scattering of photons by the water column, at the order one, given by the term L d i f shown in Figure 1a.
To consider a whole image instead of a single target pixel ρ t , the target pixel index i is now introduced in the notation, for all the concerned variables. The target pixel reflectance is decomposed as :
ρ i = j = 1 J s j a j , i = S a i ,
where a j , i is the abundance fraction of the pure material (endmember) reflectance s j in the target pixel i, and J is the number of endmembers. S is composed of the endmember spectra arranged in columns, and the vector a i holds the abundances of each endmember in pixel i.
The approximation of a limited number of influencing neighboring pixels is relevant in this investigation. For the sake of clarity, the model is developed with a neighborhood given by the N nearest pixels (usually N = 4 or N = 8 ) and could be easily generalized to a wider neighborhood. Then ρ a d j = 1 N n V t , n t ρ n , where ρ n is the reflectance of pixel n, in the neighborhood { V t , n t } of the target pixel.
Following Equation (4), the equivalent environment reflectance corresponding to a target pixel i of the sub-surface image is defined as:
ρ e n v , i = δ i j = 1 J a j , i s j + ( 1 δ i ) 1 N n V i , n i j = 1 J a j , n s j = j = 1 J n V i p n , i a j , n s j = S A p i .
The coefficient p n , i is p n , i = 1 δ i N if n V i , n i , p n , i = δ i if n = i , and else p n , i = 0 . S is the L × J matrix containing all the endmembers arranged in columns, A is the J × I abundance matrix, each of its rows corresponding to an endmember, and the I × 1 vector p i contains the coefficients p n , i .

2.5. Sub-Surface Mixing Model Considering the Adjacency Effects

The last step to obtain the sub-surface mixing model consists of expressing the sub-surface reflectance as a function of the seabed endmembers spectra and abundances. The mixing model is derived from the radiative transfer equation (Equation (1)), jointly with the development of the environment reflectance proposed in Equation (6).
The contribution of the water column to the water upward radiance, namely L u d p (see Table 1 and Figure 1a), is subtracted from Equation (2), and the equation is divided by the downward illumination at surface level E 0 in order to obtain reflectance values (here the division ⊘ is elementwise). Equation (7) is obtained as follows:
L u L u d p E 0 = E t o t π E 0 T d i r ρ t + E t o t π E 0 T d i f ρ e n v .
For each pixel i of the whole data, the following notations are introduced:
r ˜ i = [ L u L u d p ] E 0 ; k 1 , i = T d i r E t o t π E 0 ; k 2 , i = T d i f E t o t π E 0 .
By using the mixture decomposition of ρ e n v and ρ t , provided in Equation (6), the expression in (7) becomes:
r ˜ i = k 1 ( S a i ) + k 2 ( S A p i ) .
For the entire scene containing I pixels and L wavelengths, the proposed sub-surface mixing model, denoted as S M w a d j in the following, is then:
R ˜ = K 1 ( S A ) + K 2 ( S A P ) .
The L × I matrices R ˜ , K 1 , K 2 are constructed with all the spectral vectors arranged in columns, and the matrix P , which depends on the environment parameter δ i at each pixel, is an I × I matrix. The adjacency effects caused by the photons scattered by the water column are described by the additive term K 2 ( S A P ) in Equation (10).
Let us note that if the adjacency effects are not taken into consideration (i.e., assuming δ = 1 ), the matrix P becomes the identity, P = Id I × I , so the following simpler sub-surface mixing model, denoted as S M w n o a d j , is obtained:
R ˜ = ( K 1 + K 2 ) ( S A ) .

2.6. Unmixing Method

The objective of the unmixing method consists of estimating both matrices A and S . Here, the case where the water column characteristics are assumed to be known is being considered, i.e., that K 1 and K 2 are known and do not need to be determined. Such a case could be realistic when bio-optical measurements of the water column (i.e., water quality) have been carried out and when the bathymetry data are available from other means, for example using LiDAR techniques. Water composition parameters could also be estimated by using an L S inversion method that would not require the accurate knowledge of the seabed (for example S W I M code [48]), or alternatively by inversion using a fixed spectral library [36].
Then E 0 , E t o t , L u d p , T d i r , T d i f , δ can be numerically calculated for each pixel, using the O S O A A radiative transfer model, to obtain the matrices K 1 , K 2 , P . The matrix P is beforehand calculated for a chosen neighborhood size, which should depend on the pixel size and the scattering properties of the water column.
To solve the joint estimation of abundances and endmembers problem, a cost function R Q E , defined in Equation (12), is minimized.
R Q E = R ˜ o b s R ˜ F r o 2 .
The R Q E is the Frobenius norm of the difference between the observed sub-surface reflectance modified according to (8) (see also Table 1), R ˜ o b s , and the estimated one, R ˜ , given by the expression of Equation (10). R ˜ depends on the variable matrices A and S .
Since the unmixing problem is ill-posed, solutions from local minima of the cost function may be obtained, or many global minima could exist, possibly giving erroneous endmembers and abundances. Initialization is then a key step for the unmixing procedure. The application of constraints allows the range of solutions to be reduced. Generally, constraints are chosen according to some characteristics of the data which are then guaranteed to be preserved [49,50]. Two variants of N M F are developed hereafter, respectively so-called W A D J U M , corresponding to the mixing model S M w a d j expressed in Equation (10), and W U M , corresponding to the mixing model S M w n o a d j of Equation (11). The two variants are both based on the M i n i d i s c o algorithm [50], which uses an alternate projected gradient descent with an Armijo-Lin-based rule [24,51], and optional regularization, such as the sum-to-unity constraint ( S T U ) for the abundances, or the minimum dispersion constraint for the spectra [50]. The estimation of the matrices A , S is obtained following Equation (13):
[ A e s t , S e s t ] = arg min ( A , S ) R Q E ( A , S ) + Λ ( A , S ) .
Λ denotes regularization terms that can be added to the R Q E criterion. In our experiments, a soft constraint to encourage the sum-to-unity for abundances is used, Λ S T U = λ S T U ( I 1 J A I 1 I ) ( I 1 J A I 1 I ) T . The vector I 1 J is composed of the values ’one’; the regularization coefficient λ S T U is empirically set to 0.5 after various trials.
For the W A D J U M algorithm, the gradients of R Q E , relative to A and S , are detailed in Equation (14), where t denotes the transpose.
R Q E S = K 1 R ˜ o b s K 1 ( S A ) K 2 ( S A P ) A t + K 2 R ˜ o b s K 1 ( S A ) K 2 ( S A P ) ( A P ) t R Q E A = S t K 1 R ˜ o b s K 1 ( S A ) K 2 ( S A P ) + S t K 2 R ˜ o b s K 1 ( S A ) K 2 ( S A P ) P t .
In the W A D J U M unmixing algorithm, updates of A , S , and R ˜ are successively performed at each iteration. The values of S and A are restricted to vary between 0 and 1, by setting the negative values to 0 and setting to 1 the values greater than 1. The optimization is stopped after 1000 iterations or after a threshold error is attained ( 1 % of relative variation of R ˜ ).
An additive regularization term for the endmembers could also be used for example to minimize the spread of the endmembers (or the volume of the endmember simplex) as defined in [50]. However in the case of underwater unmixing, this constraint is not appropriate, due to the shape of the water column attenuation that strongly mitigates the signal between 600 nm and 700 nm (see Figure 2), resulting in possibly flat estimated spectra in this range if this constraint is applied.
Please note that the simpler variant of the underwater unmixing method for which no adjacency effects are taken into consideration, referred to as W U M , is obtained from the same algorithm, by setting P = Id I × I .

3. Analysis of Performance Based on Simulations

Simulations are carried out to evaluate the performance of the proposed W U M and W A D J U M algorithms.

3.1. Evaluation Metric

The results are evaluated using the mean spectral angle mapper (SAM) and the mean normalized R M S E , defined as follows in Equation (15):
S A M = 1 N ω ω 1 I i cos 1 ( s i t s ^ i ) ( s i t s i s ^ i t s ^ i ) ω ; N M R M S E = 1 N ω ω M M ^ F r o M F r o ω .
Here ω and i respectively refer to the realizations (different noise generation and random initialization) and to the pixels, s i and s ^ i refer to the true and estimated spectra, respectively; M and M ^ refer to the true and estimated matrix, respectively. The generic matrix M can be A , S , K 1 , or K 2 . As the solution of the unmixing method is obtained up to a permutation for the matrices A and S , the best match between true and estimated endmembers is performed before calculating the normalized R M S E .
The benchmark for evaluating the performance of the proposed algorithms is to apply the unmixing method on the seabed in the absence of a water column, in order to compare the proposed underwater unmixing algorithms with a method from the literature that considers both abundance and endmember estimation. The results of the two sub-surface unmixing algorithms W U M and W A D J U M are compared to those of the simple unmixing of the seabed with no water column, using the algorithm M i n i d i s c o [50], hereafter referred to as ` N M F , no water’, using the same random initialization. This reference allows the testing of only the influence of the perturbation due to the water column, because the gradient-based optimization used is the same as the one of W A D J U M and W U M .

3.2. Simulated Data

The influence of the water column optical properties (attenuation due to scattering and absorption by hydrosols) is accounted for using O S O A A model simulations, using the selected water turbidity and depth as inputs of O S O A A . Please note that the approximations outlined in Section 2 to establish the model given in Equation (10) were not applied on simulated data, since the radiative terms L u d p , E t o t , E 0 , T d i r , T d i f , δ were calculated using the O S O A A radiative transfer model, including all multiple reflections on the seabed and the water column. Some stationary white noise was also added to the sub-surface reflectance R ˜ to simulate environment and sensor noise, with an SNR value of 40 dB. The values of K 1 and K 2 determine the relative importance of the direct and scattered light, respectively. Figure 2 shows the spectral variation of one column of the matrices K 1 and K 2 (respectively k 1 and k 2 ), corresponding to one pixel, for various depths ranging from 1 m to 15 m, and for two water turbidities (clear water, Figure 2a; moderately turbid water, Figure 2b).
The direct attenuation vector k 1 always decreases with the depth whatever the wavelength, whereas the diffuse attenuation vector k 2 shows a more complex behavior with depth. This is because the illumination E t o t at the seabed level and the direct transmittance T d i r both systematically decrease with depth, whereas the diffuse transmittance T d i f could either increase or decrease with depth, depending on the range of seabed depth and on the turbidity of the water. Therefore, the matrix K 2 which is expressed as T d i f E t o t π E 0 , varies with depth following a compromise between depth and water turbidity. For very shallow water (e.g., seabed depth lower than 2 m), the direct attenuation term k 1 contributes more than the diffuse attenuation term k 2 , so the influence of the adjacency effect induced by neighboring pixels may not be noticeable. However, in some cases, k 2 is comparable to or even higher than k 1 , for example when the seabed depth z is greater than 10 m in clear water, and when z is greater than 2 m in the considered moderately turbid water. For those cases, the use of the model proposed in Equation (10)) for the unmixing is relevant because the adjacency effect has a significant impact on the total upward radiance.
Seabed reflectance data were generated using a random linear mixing of four endmember spectra representative of benthic species. These spectra were selected from a spectral library shown in Figure 5f (marked spectra) that had been measured during a field experiment that occurred in the Porquerolles Island, France, in September 2015. The measured spectra covered various types of seabed composition such as Caulerpa taxifolia algae, Posidonia seagrass, sand and substratum. They were sampled between 400 nm and 700 nm with a 10 nm step (31 wavelengths). The abundances were randomly generated using a maximum abundance of 85% for each endmember (no pure pixel), and some stationary white noise was added to the seabed image, in order to reach an SNR value of 40 dB at the seabed level. The hyperspectral image size is ( 100 × 24 ) pixels × 31 wavelengths .

3.3. Influence of the Consideration of the Adjacency Effects within a Seabed Unmixing Algorithm

The objective of this sub-section is to test the relevance of introducing the adjacency effects within the seabed unmixing method. To this end, a comparison is performed between the W A D J U M algorithm ( S M w a d j model), and the W U M algorithm ( S M w n o a d j model). The unmixing method is also applied in the absence of a water column, using the M i n i d i s c o algorithm [50], denoted as ` N M F , no water’ in the figures. The results of M i n i d i s c o serve as a benchmark for testing the influence of the water column on the unmixing results, using the same gradient-based optimization than the proposed NMF-based methods W A D J U M and W U M , without any influence of the water column.
The soft sum-to-unity constraint as defined in Section 2.6 is the only regularization process applied in the simulations, in order to avoid any interference with the evaluation of the influence of the adjacency effects on the unmixing algorithm.
The unmixing algorithm was initialized using biased endmembers, obtained using a random mixing of the true endmembers with other spectra from the spectral library shown in Figure 5f. Each initialization spectrum was composed of 80 % of the true endmember plus 20 % of a uniform random mixture of the true endmember and another one from the library. The abundances were then initialized using the F C L S (Fully Constrained Least-Square) algorithm [52] applied to the seabed. The algorithms were developed in the M A T L A B language. M i n i d i s c o ( N M F , no water) converged in about 4 s on a standard laptop computer (Intel(R) Core(TM) i7-461, CPU 3.00 GHz, 16 Mo RAM), whereas W A D J U M and W U M running times were longer to attain the stopping condition (about 150 s). All the algorithms were initialized in the same manner.
Two water conditions were considered, clear water ( C h l = 0.03 mg·m 3 , S M = 0.01 mg·L 1 , C D O M = 0.01 m 1 ), and moderately turbid water ( C h l = 1 mg·m 3 , S M = 1 mg·L 1 , C D O M = 0.1 m 1 ). In very turbid water conditions, due to the water attenuation, the amount of information issued from the seabed is not sufficient to envisage blind unmixing application.
Figure 3 represents the spectral angle mapper value ( S A M ), the normalized root mean square error ( N S R M S E ) on estimated endmember spectra, and the N A R M S E on estimated abundances, obtained for clear and moderately turbid waters, a pixel radius value of 0.5 m and water depth z set to [ 1 , 2 , 3 , 4 , 5 , 10 ] m. Each result is averaged over 10 realizations of noise and random initialization.
The results of the unmixing method using M i n i d i s c o , applied to the case where the water column is not present, show an angle error around 0.02 rad, a normalized R M S E around 3 % for the spectra, and around 10 % for the abundances. The presence of the water column expectedly degrades the performance of the unmixing methods. However, the range of errors remains acceptable for both algorithms; as an example, in the clear water case the angle error is around 0.03 rad and the normalized R M S E is around 6 % for spectra, and the N A R M S E is around 12 % for abundances. In clear water, the difference between W U M and W A D J U M is not tremendous, because the direct term K 1 is greater than the diffuse term K 2 , so the adjacency effect has small relative influence. In the case of moderately turbid water, the W U M algorithm (no adjacency) gives very poor performance from z higher than 4 m, thus making this unmixing method not applicable (relative errors greater than 30 % for spectra and greater than 40 % for abundances), whereas the performance of W A D J U M remains about the same as the one obtained in clear water context. Standard deviations are comparable for the W U M and W A D J U M algorithms in the case of clear water, whereas the R M S E obtained using the W U M algorithm appears to be unstable for moderately turbid waters at the intermediate seabed depth of 5 m. Therefore, the proposed W A D J U M algorithm allows efficient performance of the unmixing of the seabed with specifications close to those obtained when the unmixing is carried out in the absence of a water column. In moderately turbid water, it appears necessary that the unmixing method takes into consideration the adjacency effects to estimate the seabed endmembers and corresponding abundances.

3.4. Analysis of the Robustness of the Unmixing Methods with Respect to a Lack of Knowledge of the Water Column Properties

In the previous sections, K 1 and K 2 were supposed to be perfectly known. However, in real-world conditions, knowledge of the water attenuation is not perfect.
It is reminded here that the values of K 1 , K 2 and δ (thus the matrix P ) depend both on the water quality, which is more or less homogeneous in small areas, and on the bathymetry, which could potentially vary from one pixel to another. To examine the impact of an error made on the water column optical properties, the true values of K 1 , K 2 and δ are replaced by biased ones in the inputs of W U M and W A D J U M algorithms. The way in which the errors on K 1 , K 2 and δ are generated has a significant influence on the results. Although it is not possible to test all configurations, two cases are analyzed hereafter.

3.4.1. Error on the Bathymetry

Here the performance of the algorithms is evaluated when a random error on the bathymetry is applied, for two values of the true depth, z t r u e = 5 . 5 m and z t r u e = 9 . 5 m, in the case of moderately turbid water. The biased values of the bathymetry are obtained using a uniform random generation of the seabed depth, U [ 0 . 5 , 0 . 5 ] m, centered on the true values. The corresponding biased values of K 1 , K 2 and δ , used for the input of the algorithms, are obtained by interpolating between the values given by O S O A A simulations for z = 5 m and z = 6 m (resp. z = 9 m and z = 10 m). The performance of the unmixing methods is presented in Table 3.
When z t r u e = 5 . 5 m, the relative error on K 1 and K 2 due to the error on bathymetry is not very high ( N K 1 R M S E = 6 % , N K 2 R M S E = 5 % ). However the unmixing results remain satisfactory for W A D J U M only ( N S R M S E = 4 % , N A R M S E = 14 % ), whereas the W U M algorithm provides unsatisfactory estimate for spectra ( N S R M S E = 22 % ) and for abundances ( N A R M S E = 30 % ).
For the deepest seabed case ( z t r u e = 9 . 5 m), the random error on the depth leads to a significant error on the attenuation matrices, N K 1 R M S E = 47 % and N K 2 R M S E = 18 % . This high error on the water attenuation results in degraded performance for both W U M ( N S R M S E = 41 % , N A R M S E = 60 % ) and W A D J U M ( N S R M S E = 22 % , N A R M S E = 29 % ), although the estimation remains better for W A D J U M .
In any case, the W A D J U M algorithm still provides better results than the W U M algorithm, so the relevance of taking in consideration the adjacency effects for the unmixing is confirmed. In moderately turbid water, the W A D J U M algorithm is robust to a uniform error of amplitude 1 m on the bathymetry at a depth of 5 . 5 m, but this error has a greater impact at a depth of 9 . 5 m, thus limiting the field of application of the unmixing method in the case of biased knowledge of the bathymetry.

3.4.2. Error on the Knowledge of the Water Quality

The objectives of this paragraph are (i) to study the consequences of a bias on the knowledge of water quality, (ii) to verify if scenes with constant bathymetry are easier to map than scene with varying bathymetry.
The test image contains five areas of each 100 × 24 pixels and 31 spectral bands, each area corresponding to one fixed depth, respectively z = 1 m , 2 m , 3 m , 4 m , 5 m . The results obtained by the W U M and W A D J U M algorithms, applied to each constant depth area separately, are compared to those obtained when the unmixing method is applied to the entire image composed of all five areas. This makes it possible to evaluate the method both when the depth is not constant and when the depth is fixed. For the entire image, the estimation of abundances is obtained globally but evaluated in each area separately, whereas the estimation and evaluation of endmember spectra is only global. The true water quality is the moderately turbid one. The erroneous values of K 1 , K 2 , δ are generated by interpolating between their respective values obtained for clear and moderately turbid water with a fixed proportion (resp. 5 % and 95 % ), for each depth. The relative R M S E errors between true and biased values of K 1 , K 2 are presented in Table 4.
It can be seen that the normalized error on the direct term, N K 1 R M S E , becomes significant when the depth increases, while N K 2 R M S E does not vary much with the depth. The results are reported in Figure 4.
For the case of separate constant seabed depth areas, due to the increasing error on K 1 in deepest parts, the error on the endmember spectra estimate, N S R S M E , also increases with depth, although it is lower using W A D J U M than W U M . The endmember spectra cannot be accurately retrieved due to the bias on water attenuation. The abundance estimate on the separate areas remains quite good for all the depths despite the distorted retrieval of the endmember spectra, as observed in Figure 4. This could be explained by the constant depth in one area, which results in the same distortion for all spectra, thus still corresponding to distinguishable classes of benthic species.
On the contrary, for the entire image, the error on endmember spectra estimate, N S R S M E , shows an intermediate value between the performance obtained at lower and higher bathymetries on the separate parts. To interpret this result, it should be reminded that the parts of the global image corresponding to the deepest areas have the smallest reflectance values, due to the attenuation induced by the water column, so these parts contribute more weakly to the R Q E , which is the criterion that is optimized during the unmixing process. Therefore, the global optimization essentially involves the shallowest parts of the entire image. Thus, the global endmember estimate remains relatively good. For the entire image, the errors on abundances obtained using W U M and W A D J U M both increase slightly on the shallowest part when using the global unmixing method; this could be due to the degradation of the global endmember estimates at z = 1 m. W A D J U M is more performing than W U M in the deepest parts.
In conclusion, a small variation in the water quality knowledge may lead to quite a large bias on K 1 at the largest seabed depths. The endmember estimation degrades with this bias and is not reliable for a relative error of the water attenuation greater than 14 % . The abundance estimation is more robust, the error remains acceptable (maximum 13 % N A R M S E ) when using W A D J U M in the case of biased water attenuation knowledge. The difference between W A D J U M and W U M performance may be not significant for constant depth areas. The performance of the abundance estimates could deteriorate when observing scenes characterized by highly variable seabed depths.

4. Application to Real Data

An airborne data campaign has been conducted by the Hytech Imaging company [53] on 13 September 2017, using the H y S p e x sensor [54], over the Porquerolles Island coast (France). The hyperspectral images were corrected to obtain georeferenced reflectance images. The pixel size of the observed images is 0.5 m × 0.5 m, the wavelengths have been resampled to achieve a 10 nm resolution, and the images were analyzed in the spectral range from 410 nm to 700 nm. Sun glint [47] and the air/water interface [22] were corrected to finally obtain the sub-surface remote-sensing reflectance image R ˜ (in sr 1 ) by subtracting the water upward reflectance. In situ measurements were performed with the help of the PNPC (Parc National de Port Cros) [55] at six stations in the Alicaster Bay, where the seabed depth ranges from 2 m to 20 m. In situ data consisted of measuring chlorophyll-a ( C h l ), suspended matter ( S M ) and colored dissolved organic matter ( C D O M ) concentrations. A diver brought back samples of the benthic species from the seabed to measure their spectral signatures (Figure 5d,f). A Remotely Operated underwater Vehicle (ROV) data campaign has been carried out by the French Institute for Marine Exploration (Ifremer) [56] over a transect of 1500 × 2 m surface, with a high spatial resolution ( 2 . 3 cm), to acquire R G B images at the seabed level [57]. The images have been classified in three classes, namely sand, Posidonia, and Caulerpa taxifolia. The supervised classification method is based on color information, using a weighted combination of Euclidian and spectral angle distance, and reference spectra representative of each class. This classification allows the generation of an abundance map of 0 . 5 m spatial resolution, which serves as a ground truth data for the evaluation of the unmixing results. A colored composition of the sub-surface image, with the vehicle’s path, is shown in Figure 5a, together with the bathymetry and an example of R O V underwater image containing sand, Posidonia and Caulerpa taxifolia (Figure 5b,c).
The mean measured values of the concentrations of the water constituents in the Alicaster Bay were C h l = 0.03 mg·m 3 , S M = 0.5 mg·L 1 , C D O M = 0.05 m 1 . The bathymetry was provided by the L I T T O 3 D database, available in [58], supported by the I G N (French National Institute for Geographical and Forestry Information) and the S H O M (Marine Hydrographic and Oceanographic Service). Please note that the entire French coast bathymetry is available with a precision of 95 % , and 1 m spatial resolution. The values of the water attenuation K 1 and K 2 , and δ , have been calculated beforehand using the OSOAA radiative transfer model, the mean measured values of the concentrations of water constituents, C h l , S M , C D O M , and the bathymetry.

4.1. Initialization of the Unmixing

Initial matrices S i n i t and A i n i t must be determined for the unmixing procedure. In order to find such matrices, two steps are performed.
Step one: a first sea column correction is carried out from a classical L S method (lsqcurvefit function of M A T L A B ), using a chosen library S 0 , and the model S M w n o a d j of Equation (11), to obtain the coefficient matrix C e s t . Although the spectra in S 0 should be representative of the local benthic habitat, they are not supposed to be the exact endmembers. Consequently, the sum-to-unity constraint is ignored in the inversion, in order to permit scale factors for the spectra. A preliminary estimate of the seabed reflectance, given by X i n v = S 0 * C e s t , is then recovered.
Step two : the initial matrix S i n i t is obtained by applying the V C A algorithm [59] to X i n v , and the initial abundance matrix A i n i t is obtained using the F C L S algorithm [52] and S i n i t .

4.2. Experiments on the Validation Zone

The validation data have been divided in two parts, Part 1 and Part 2, as shown in Figure 5a, in order to avoid the use of big matrices for calculation. A spectral library has been acquired in the Porquerolles area in June 2016, September 2016 and September 2017. Different spectra representative of sand, substratum and Posidonia are shown in Figure 6.
It can be seen that the Posidonia class can be represented by various spectral signatures, as well as sand or substratum. Moreover, organic surface deposits such as waste stems can be found on dark sand and substratum, so their reflectances are very close to those of brown Posidonia or shell colonized Posidonia. Consequently, in large areas such as the one covered by the robot transect, some confusion can be expected between Posidonia and dark sand or substratum classes. In the following, only estimated abundances ( A e s t ) and coefficients ( C e s t ) are shown, the estimated spectra being distorted due to the approximative knowledge of the water attenuation.

4.2.1. Abundance Estimation Performed on Part 1

In this area, only Posidonia seagrass and sand are present. We consider two reflectance spectra (green and brown) for the Posidonia class, and the clear sand spectrum for the sand, which is locally light and homogeneous (see Figure 6). The ground truth abundances A G T , the L S estimated coefficients C e s t , and the W A D J U M estimated abundances A e s t are respectively shown in Figure 7a–c, as a function of the pixel index. The abundances/coefficients of the two P o s i d o n i a classes have been fused. The N A R M S E errors are calculated for the three unmixing methods: V C A + F C L S performed on the reconstructed seabed X i n v with no water column, and W U M / W A D J U M performed on the sub-surface reflectance image R ˜ .
It can be seen that for both L S estimation and unmixing methods, some sand is detected where the ground truth indicates only Posidonia. This is due to the variability of the Posidonia reflectance, which can be colonized by shells, so with a spectral signature close to the sand’s one. The distributions of the abundances of the two classes estimated using W A D J U M better fit the ground truth than the L S coefficients, which shows the interest of the unmixing method. The abundance errors obtained for the three unmixing methods are N A R M S E = 22 % using V C A + F C L S performed on X i n v , N A R M S E = 19 % using W U M or W A D J U M performed on R ˜ . As two classes are present, the error is the same for the sand and for the P o s i d o n i a . The error is not calculated for the coefficients C e s t , because they are not normalized.

4.2.2. Abundance Estimation Performed on Part 2

In the Part 2 of the validation area there is a C a u l e r p a colonization surface, located in craters from ancient projectile testing. The bathymetry varies from 6 m to 14 m, and is around 10 m in the C a u l e r p a location. As the monitoring of the C a u l e r p a is an important issue for the PNPC, we focus on the estimation of such a benthic habitat. P o s i d o n i a and sand/substratum are in the same class. Results are presented in Figure 8. Please note that the ground truth shows incomplete classification at pixel ca. 85 (Figure 8a).
The distribution of C a u l e r p a given by the L S coefficients C e s t is more dispersed than the distribution given by A e s t , obtained using W A D J U M . The difference between the initialization using V C A + F C L S performed on X i n v , and the output of W U M and W A D J U M performed on R ˜ , is not tremendous for the estimation of C a u l e r p a abundance: N A R M S E = 15 % using V C A + F C L S , N A R M S E = 13 % using W U M , N A R M S E = 12 % using W A D J U M .
In conclusion, the validation shows a real improvement using an unmixing strategy instead of L S coefficient estimation, and a slight improvement using W U M or W A D J U M performed on the sub-surface images compared to the use of V C A + F C L S performed on the reconstructed seabed with no water column. The difference between W U M and W A D J U M performances is very low, according to the simulations results in clear water condition such as the clear water in Porquerolles area (see Section 3.3). Let note that both V C A + F C L S and W U M results should degrade in moderately turbid conditions, because these two methods do not take into account the adjacency effects. Furthermore, the water quality is not accurately known, further reducing the difference in performance between the W U M and W A D J U M unmixing methods (see Section 3.4.2).

4.3. Experiment on Another Area

A small area of 200 × 530 pixels (100 m × 265 m), delimited by a red rectangle in the Figure 5a, has been selected as shown in Figure 9b. This area shows a diversified seabed, and the bathymetry is in the range between 2 m and 6 m (Figure 9d). The endmembers used for the L S inversion were dark green P o s i d o n i a , sand, substratum and photophilic algae, given in Figure 5f. Photophilic algae consist of a rich variety of algae conditioned by light penetration, often located on rocky areas. They constitute an important environmental issue. The environment parameter δ , K 1 , K 2 values and the water reflectance L u d p E 0 were calculated in the same way as in Section 4.2. Because the size of the matrix P is too large for computation on M A T L A B , the image has been split into five contiguous sub-images to run the W A D J U M algorithm, and the obtained abundances have been reassembled (Figure 9c). This presupposes the hypothesis of the same classes of pure materials being made for the five sub-images, and also that some variations in the spectra representative of one class could be permitted. Please note that the imperfect knowledge of the water column (i.e., the water quality) resulted in distorted endmember spectra as mentioned in Section 3 (not shown). As the substratum and brown P o s i d o n i a Classes may overlap (see Section 4.2), a single class, denoted as Subs/Pos in Figure 9a,c is associated with both brown Posidonia and substratum classes.
The abundances estimated using W A D J U M , as shown in Figure 9c, can be qualitatively evaluated by comparison with the sub-surface image (Figure 9b). The sand is correctly mapped. The green P o s i d o n i a seagrass is observed in the middle of the area identified as substratum/brown P o s i d o n i a , and photophilic algae are located on the edge of the area of settlement of the P o s i d o n i a and substratum. In conclusion, the qualitative examination of the seabed mapping indicates a distribution of the benthic habitat that is consistent with the observed sub-surface image.
The comparison between the coefficients given by C e s t (Figure 9a) and the abundances obtained in A e s t (Figure 9c), shows some differences in the spatial distribution of benthic habitat (apart from the scale, which is not the same for abundances and coefficients), circled in red. The green P o s i d o n i a is erroneously detected on the seabed right corner of the coefficients in C e s t , and the substratum/brown P o s i d o n i a could be confused with sand in C e s t .

5. Discussion

The first question addressed in this work is: is it relevant to perform seabed unmixing, due to attenuation by the water column?
An accurate model of the radiative transfer was necessary to overcome the low magnitude of the signal issued from the seabed at the sub-surface level. The O S O A A model allowed us to select the radiative components from the complete equations of the radiative transfer inside the water column that were sufficiently impacting the sub-surface signal to be distinguished from noise. Although many non-linear interactions such as multiple reflections are present in the water column, an analysis of those contribution to the total bottom illumination showed that those non-linear phenomena do not significantly influence the upward radiance at the sub-surface level, for depths comprised between 1 m and 10 m. In consequence the only considered interactions of the bottom with the water column taken into consideration in this study are the adjacency effects. It should be highlighted that the consideration of adjacency effects within the unmixing method is the strength and the main originality of the proposed algorithm, relatively to previous studies.
Simulation results show that the W A D J U M unmixing algorithm (taking into account adjacency effects) is able to achieve satisfactory results, which are close to those obtained by unmixing the seabed without a water column, when the water column properties are perfectly known. However, the results depend on the abundance and endmember initialization matrices. A method for determining appropriate initialization matrices for real data has been proposed in Section 4.
The second question addressed is: when should adjacency effects be taken into consideration in the unmixing process?
The application of the unmixing technique using W U M (which ignores the adjacency effects), is fairly inefficient when the bathymetry is greater than 4 m for moderately turbid water (spectra estimation degrades from 2 m), while in clear water the unmixing performed using W U M and W A D J U M provided a similar performance. Another important parameter is the spatial resolution of the image, which depends on the sensor. The environment parameter δ , which gives the relative importance of adjacency effects, significantly depends on the bathymetry and on the size of the pixel. For a pixel radius larger than 1 m and shallow moderately turbid waters ( z 1 m), the environment parameter is δ 1 , meaning that the photons contributing to the upward radiance mainly originate from the seabed target itself. At a 5 m bathymetry, δ 0 . 75 , for a pixel radius of 1 m, and δ 0 . 65 for a pixel radius of 0 . 5 m. From these considerations and from the results shown in Figure 3, we conclude that it is necessary to take into account the adjacency effects for moderately turbid waters, bathymetry higher than 2 m, and for spatial image resolution lower than 2 m. It should also be mentioned that due to the complexity of the algorithm, the unmixing technique is not applicable at once on big areas. In our work we considered at most 2 × 10 4 pixels to process with the matrix P .
In the case of biased knowledge of bathymetry and water quality, the performances of the W A D J U M and W U M algorithms both decrease, and become more similar. The endmember estimation is corrupted by the erroneous water attenuation and is not reliable; however, the abundance estimation appears to be robust and remains acceptable using W A D J U M ( 13 % of relative error in our simulations), especially for constant depth areas, and for the water attenuation relative error in the range of 14 % or less.
Both unmixing algorithms were applied to real airborne hyperspectral data from the Porquerolles Island coast (France). The estimated abundance maps were compared to a ground truth given by color images issued from a R O V , in a depth range of 5.5 m–14 m. For these real data, the clear water condition induced no significant difference between W U M and W A D J U M performances; this is in accordance with the results obtained in the simulation using synthetic data, which showed that in clear water, the performances of W U M and W A D J U M were similar up to 10 m depth, for abundance estimation. Furthermore, the imperfect knowledge of the water quality also reduces the difference between W U M and W A D J U M performances. On the other hand, the comparison of the estimated abundances A e s t with the coefficients C e s t obtained using a fixed spectral library and a least-square estimation, showed the interest of an unmixing method to map the benthic habitat, in comparison with L S coefficient estimation using a fixed spectral library. In the validation area, the comparison with initialization matrices showed a small improvement using W A D J U M , which could be due to imperfect knowledge of the water bio-optical properties, and also to the clear water conditions.

6. Conclusions and Perspectives

In conclusion, we propose a strategy able to perform seabed unmixing for coastal areas analysis, with depths ranging from 1 m to 10 m, for clear or moderately turbid water, as long as enough knowledge of the water bio-optical properties is available. For clear waters or for pixel sizes higher than 2 m, the W U M algorithm can be used, while for moderately turbid waters and pixel sizes lower than 2 m, W A D J U M , which takes into consideration the adjacency effects in the water column, should be used. The interest of such a method resides in performing a precise seabed analysis in relatively small areas, while for very large areas coefficient estimation is more adapted, due to the calculus time.
Future works could be dedicated to improving the seabed unmixing method. First, the running time could be reduced by exploiting the sparsity of the matrix P that contains the environment parameter. Secondly, a method to improve the optimization strategy could be searched. For example, the Alternating Direction Method of Multipliers ( A D M M ) strategy has shown to be able to improve unmixing in difficult cases [43]. Such a strategy could be explored in the future. To improve the seabed unmixing method, the model could also evolve. Based on the presented analysis, the more adequate direction seems to take into consideration the spectral variability, as benthic habitat spectral signatures are subjected to seasonal changes.
Since the accuracy of the water composition knowledge is important for the quality of the unmixing, all methods that could improve this knowledge would be relevant to enhance the performance of W A D J U M . Last, efforts for constructing a parametric model for the water attenuation could be made in the context of the proposed radiative transfer model with adjacency effects taken into consideration.

Author Contributions

Conceptualization and methodology, M.G. and Y.D.; software, M.G. and L.J.; validation, M.G., Y.D., M.C., A.M.; formal analysis, M.G.; investigation, all authors; resources A.M., M.G., M.C., X.L. and B.L.; data curation, A.M., X.L. and B.L.; writing—original draft preparation, visualization, M.G.; supervision, M.G. and V.S.; project administration, M.G.; funding acquisition, V.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the French Defense Agency, Direction Générale de l’Armement ( D G A ), as an A S T R I D / A N R project supervised by the French national research agency Agence Nationale de la Recherche, grant HypFoM ANR-15ASTR-019. The APC was funded by IRAP-CNRS (HypFom Project).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hochberg, E.J. Remote Sensing of Coral Reef Processes In Coral Reefs: An Ecosystem in Transition; Dubinsky, Z., Stambler, N., Eds.; Springer: Dordrecht, The Netherlands, 2011; pp. 25–35. [Google Scholar]
  2. Chauvaud, S.; Bouchon, C.; Maniere, R. Remote sensing techniques adapted to high resolution mapping of tropical coastal marine ecosystems (coral reefs, seagrass beds and mangrove). Int. J. Remote Sens. 1998, 19, 3625–3639. [Google Scholar] [CrossRef]
  3. Jaubert, J.; Chisholm, J.; Minghelli-Roman, A.; Marchioretti, M.; Morrow, J.; Ripley, A. Re-evaluation of the extent of caulerpa taxifolia development in the northern mediterranean using airborne spectrographic sensing. Mar. Ecol. Prog. Ser. 2003, 263, 75–82. [Google Scholar] [CrossRef] [Green Version]
  4. Lee, Z.; Carder, K. Effect of spectral band numbers on the retrieval of water column and bottom properties from ocean color data. Appl. Opt. 2002, 41, 2191–2201. [Google Scholar] [CrossRef] [PubMed]
  5. Minghelli-Roman, A.; Chisholm, J.; Marchioretti, M.; Jaubert, J. Discrimination of coral reflectance spectra in the red sea. Coral Reef 2002, 21, 307–314. [Google Scholar] [CrossRef]
  6. Lee, Z.; Weidemann, A.; Arnone, R. Combined effect of reduced band number and increased bandwidth on shallow water remote sensing: The case of worldview 2. IEEE Trans. Geosci. Remote Sens. 2013, 51, 2577–2586. [Google Scholar] [CrossRef]
  7. Hochberg, E.; Atkinson, M.; Andréfouët, S. Spectral reflectance of coral reef bottom-types worldwide and implications for coral reef remote sensing. Remote Sens. Environ. 2003, 85, 159–173. [Google Scholar] [CrossRef]
  8. Botha, E.; Brando, V.; Anstee, J.; Dekker, A.; Sagar, S. Increased spectral resolution enhances coral detection under varying water conditions. Remote Sens. Environ. 2013, 131, 247–261. [Google Scholar] [CrossRef]
  9. Hochberg, E.; Atkinson, M. Capabilities of remote sensors to classify coral, algae, and sand as pure and mixed spectra. Remote Sens. Environ. 2003, 85, 174–189. [Google Scholar] [CrossRef]
  10. Hedley, J.; Roelfsema, C.; Phinn, S.; Mumby, P. Environmental and sensor limitations in optical remote sensing of coral reefs: Implications for monitoring and sensor design. Remote Sens. 2012, 4, 271–302. [Google Scholar] [CrossRef] [Green Version]
  11. Hedley, J.; Roelfsema, C.; Phinn, S. Efficient radiative transfer model inversion for remote sensing applications. Remote Sens. Environ. 2009, 113, 2527–2532. [Google Scholar] [CrossRef]
  12. Lee, Z.; Carder, K.; Mobley, C.; Steward, R.; Patch, J. Hyperspectral remote sensing for shallow waters. I. A semianalytical model. Appl. Opt. 1998, 37, 6329–6338. [Google Scholar] [CrossRef] [PubMed]
  13. Maritorena, S.; Morel, A.; Gentili, B. Diffuse reflectance of oceanic shallow waters - influence of water depth and bottom albedo. Limnol. Oceanogr. 1994, 39, 1689–1703. [Google Scholar] [CrossRef]
  14. Mobley, C. Light and Water: Radiative Transfer in Natural Waters; Academic Press: Cambridge, MA, USA, 1994. [Google Scholar]
  15. Chami, M.; Lafrance, B.; Fougnie, B.; Chowdhary, J.; Harmel, T.; Waquet, F. Osoaa: A vector radiative transfer model of coupled atmosphere-ocean system for a rough sea surface application to the estimates of the directional variations of the water leaving reflectance to better process multi-angular satellite sensors data over the ocean. Opt. Express 2015, 23, 27829. [Google Scholar] [PubMed] [Green Version]
  16. Dekker, A.G.; Phinn, S.R.; Anstee, J.; Bissett, P.; Brando, V.E.; Casey, B. Intercomparison of shallow water bathymetry, hydro-optics, and benthos mapping techniques in australian and caribbean coastal environments. Limnol. Oceanol. Methods 2011, 9, 396. [Google Scholar] [CrossRef] [Green Version]
  17. Jay, S.; Guillaume, M. A novel maximum likelihood based method for mapping depth and water quality from hyperspectral remote-sensing data. Remote Sens. Environ. 2014, 147, 121–132. [Google Scholar] [CrossRef]
  18. Jay, S.; Guillaume, M.; Minghelli, A.; Deville, Y.; Chami, M.; Lafrance, B.; Serfaty, A. Hyperspectral remote sensing of shallow waters: Considering environmental noise and bottom intra-class variability for modeling and inversion of water reflectance. Remote Sens. Environ. 2017, 200, 352–367. [Google Scholar] [CrossRef]
  19. Jay, S.; Guillaume, M. Regularized estimation of bathymetry and water quality using hyperspectral remote sensing. Int. J. Remote Sens. 2016, 37, 263–289. [Google Scholar] [CrossRef]
  20. Thompson, D.; Hochberg, E.; Asner, G.; Green, R.; Knapp, D.; Gao, B.-C.; Garcia, R.; Gierach, M.; Lee, Z.; Maritorena, S.; et al. Airborne mapping of benthic reflectance spectra with bayesian linear mixtures. Remote Sens. Environ. 2017, 200, 18–30. [Google Scholar] [CrossRef]
  21. Goodman, J.A. Hyperspectral Remote Sensing of Coral Reefs: Deriving Bathymetry, Aquatic Optical Properties and a Benthic Spectral Unmixing Classification Using Aviris Data in the Hawaiian Islands. Ph.D. Dissertation, Hydrologic Sciences, Department of Land, Air and Water Resources, University of California, Davis, CA, USA, 2004. [Google Scholar]
  22. Lee, D.D.; Seung, H.S. Learning the parts of objects by nonnegative matrix factorization. Nature 1999, 401, 788–791. [Google Scholar] [CrossRef]
  23. Lee, D.D.; Seung, H.S. Algorithms for non-negative matrix factorization. Adv. Neural Inform. Process. Syst. 2001, 13, 556–562. [Google Scholar]
  24. Cichocki, A.; Zdunek, R.; Phan, A.H.; Amari, S.-I. Nonnegative Matrix and Tensor Factorizations; Wiley: Chichester, UK, 2009. [Google Scholar]
  25. Bioucas-Dias, J.M.; Plaza, A.; Dobigeon, N.; Parente, M.; Du, Q.; Gader, P.; Chanussot, J. Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches. IEEE J. Sel. Top. Appl. Earth Observat. Remote Sens. 2012, 5, 354–379. [Google Scholar] [CrossRef] [Green Version]
  26. Heylen, R.; Parente, M.; Gader, P. A review of nonlinear hyperspectral unmixing methods. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2014, 7, 1844–1868. [Google Scholar] [CrossRef]
  27. Meganem, I.; Deville, Y.; Hosseini, S.; Deliot, P.; Briottet, X. Linear quadratic blind source separation unsing nmf to unmix urban hyperspectral images. IEEE Trans. Geosci. Remote Sens. 2014, 62, 1822–1833. [Google Scholar]
  28. Eches, O.; Dobigeon, N.; Tourneret, J.-Y. Enhancing hyperspectral image unmixing with spatial correlations. IEEE Trans. Geosci. Remote Sens. 2011, 49, 4239–4247. [Google Scholar] [CrossRef] [Green Version]
  29. Heylen, R.; Scheunders, P. A multilinear mixing model for nonlinear spectral unmixing. IEEE Trans. Geosci. Remote Sens. 2015, 54, 240–251. [Google Scholar] [CrossRef]
  30. Revel, C.; Deville, Y.; Achard, V.; Briottet, X. Inertia-constrained pixel-by-pixel nonnegative matrix factorisation: A hyperspectral unmixing method dealing with intra-class variability. Remote Sens. 2017, 10, 1706. [Google Scholar] [CrossRef] [Green Version]
  31. Thouvenin, P.-A.; Dobigeon, N.; Tourneret, J.-Y. Hyperspectral unmixing with spectral variability using a perturbed linear mixing model. IEEE Trans. Signal Process. 2016, 64, 525–538. [Google Scholar] [CrossRef] [Green Version]
  32. Halimi, A.; Honeine, P.; Bioucas-Dias, J.-M. Hyperspectral unmixing in presence of endmember variability, nonlinearity, or mismodeling effects. IEEE Trans. Image Process. 2016, 25, 4565–4579. [Google Scholar] [CrossRef] [Green Version]
  33. Goodman, J.; Ustin, S. Classification of benthic composition in a coral reef environment using spectral unmixing. J. Appl. Remote Sens. 2007, 1, 011501. [Google Scholar]
  34. Lee, Z.; Carder, K.; Mobley, C.; Steward, R.; Patch, J. Hyperspectral remote sensing for shallow waters. 2. Deriving bottom depths and water properties by optimization. Appl. Opt. 1999, 38, 3831–3843. [Google Scholar] [CrossRef] [Green Version]
  35. Torres-Madronero, M.; Velez-Reyes, M.; Goodman, A. Underwater unmixing and water optical properties retrieval using hyciat. In Proceedings of SPIE: Imaging Spectrometry XIV; SPIE Optical Engineering + Applications: San Diego, CA, USA, 2009; Volume 7457. [Google Scholar]
  36. Marcello, J.; Eugenio, F.; Martín, J.; Marqués, F. Seabed Mapping in Coastal Shallow Waters Using High Resolution Multispectral and Hyperspectral Imagery. Remote Sens. 2018, 10, 1208. [Google Scholar] [CrossRef] [Green Version]
  37. Guillaume, M.; Michels, Y.; Jay, S. Joint estimation of water column parameters and seabed reflectance combining ML and unmixing algorithm. In Proceedings of the 2015 7th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), Tokyo, Japan, 2–5 June 2015. [Google Scholar]
  38. Krupa, A.J.D.; Samiappan, D.; Hemalatha, V. Techniques for seabed mapping using underwater hyperspectral imaging: A survey. Int. J. Pure Appl. Math. 2018, 118, 11–30. [Google Scholar]
  39. Feng, L.; Hu, C. Cloud adjacency effects on top-of-atmosphere radiance and ocean color data products: A statistical assessment. Remote Sens. Environ. 2016, 174, 301–313. [Google Scholar] [CrossRef]
  40. Bulgarelli, B.; Kiselev, V.; Zibordi, G. Simulation and analysis of adjacency effects in coastal waters: A case study. Appl. Opt. 2014, 53, 1523–1545. [Google Scholar] [CrossRef] [PubMed]
  41. Sei, A. Analysis of adjacency effects for two Lambertian half-spaces. Int. J. Remote Sens. 2007, 28, 1873–1890. [Google Scholar] [CrossRef]
  42. Santer and Schmechtig. Adjacency effects on water surfaces: Primary scattering approximation and sensitivity study. Appl. Opt. 2000, 39, 361–375. [Google Scholar] [CrossRef] [PubMed]
  43. Wang, X.; Zhong, Y.; Zhang, L.; Xu, Y. Blind Hyperspectral Unmixing Considering the Adjacency Effect. IEEE Trans. Geosci. Remote Sens. 2019, 57, 6633–6649. [Google Scholar] [CrossRef]
  44. Burazerović, D.; Heylen, R.; Geens, B.; Sterckx, S.; Scheunders, P. Detecting the adjacency effect in hyperspectral imagery with spectral unmixing techniques. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2013, 6, 1070–1078. [Google Scholar] [CrossRef]
  45. Chami, M.; Lenot, X.; Guillaume, M.; Lafrance, B.; Briottet, X.; Minghelli, A.; Jay, S.; Deville, Y.; Serfaty, V. Analysis and quantification of seabed adjacency effects in the sub-surface upward radiance in shallow waters. Opt. Express 2019, 27, A319–A338. [Google Scholar] [CrossRef]
  46. Guillaume, M.; Juste, L.; Lenot, X.; Deville, Y.; Lafrance, B.; Chami, M.; Jay, S.; Minghelli, A.; Briottet, X.; Serfaty, V. NMF hyperspectral unmixing of the sea bottom: Influence of the adjacency effects, model and method. In Proceedings of the 2018 9th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), Amsterdam, The Netherlands, 23–26 September 2018. [Google Scholar]
  47. Hedley, H.A.; D, J.; Mumby, P. Technical note: Simple and robust removal of sun glint for mapping shallow-water benthos. Int. J. Remote Sens. 2005, 26, 2107–2112. [Google Scholar] [CrossRef]
  48. Sicot, G.; Lennon, M.; Corman, D.; Gauthiez, F. Estimation of the sea bottom spectral reflectance in shallow water with hyperspectral data. In 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS); IEEE: Milan, Italy, 2015; pp. 2311–2314. [Google Scholar] [CrossRef]
  49. Deville, Y. From separability/identifiability properties of bilinear and linear-quadratic mixture matrix factorization to factorization algorithms. Digit. Signal Process. 2019, 87, 21–33. [Google Scholar] [CrossRef]
  50. Huck, A.; Guillaume, M.; Blanc-Talon, J. Minimum dispersion constrained nonnegative matrix factorization to unmix hyperspectral data. IEEE Trans. Geosci. Remote Sens. 2010, 48, 2590–2600. [Google Scholar] [CrossRef]
  51. Lin, C.-J. Projected gradient methods for non-negative matrix factorization. Neural Comput. 2007, 19, 2756–2779. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Heinz, D.; Chang, C.-I. Fully constrained least-squares linear spectral mixture analysis method for material quantification in hyperspectral imagery. IEEE Trans. Geosci. Remote Sens. 2001, 29, 529–545. [Google Scholar] [CrossRef] [Green Version]
  53. Hytech Imaging. Available online: https://hytech-imaging.fr/ (accessed on 4 July 2017).
  54. Hyspex. Available online: https://www.hyspex.com/hyspex-products/ (accessed on 10 December 2016).
  55. Parc National de Port Cros. Available online: http://www.portcros-parcnational.fr (accessed on 10 December 2016).
  56. Institut français de recherche pour l’exploitation de la mer. Available online: https://wwz.ifremer.fr (accessed on 4 July 2017).
  57. Rigaud, V.; le Rest, E.; Marce, L.; Maniere, E.C.; Simon, D.; Peuch, A.; Perrier, M. VORTEX: Versatile and open subsea robot for technical experiment: Prototyping software architecture for the next AUV and ROV generation. In Proceedings of the 4th International Offshore and Polar Engineering Conference, Osaka, Japan, 10–15 April 1994. [Google Scholar]
  58. LITTO3D. Available online: https://www.geoportail.gouv.fr/donnees/litto3d (accessed on 15 December 2017).
  59. Nascimento, J.M.; Bioucas-Dias, J.M. Vertex component analysis: A fast algorithm to unmix hyperspectral data. IEEE Trans. Geosci. Remote Sens. 2005, 43, 898–910. [Google Scholar] [CrossRef] [Green Version]
Figure 1. (a) Various components of the light that contributes to the sub-surface radiance (from [45]). The reader is referred to Table 1 for the definition of the notations used in the Figure. (b) Environment vector δ , which represents the weight of the target pixel contributing to the diffuse radiance, for a target radius of 0 . 5 m, for various seabed depths ranging from 1 m to 10 m, and moderately turbid water conditions, corresponding to a chlorophyll-a concentration C h l = 1 mg m 3 , a suspended matter concentration S M = 1 mg L 1 , and a color dissolved organic matter concentration C D O M = 0.1 m 1 .
Figure 1. (a) Various components of the light that contributes to the sub-surface radiance (from [45]). The reader is referred to Table 1 for the definition of the notations used in the Figure. (b) Environment vector δ , which represents the weight of the target pixel contributing to the diffuse radiance, for a target radius of 0 . 5 m, for various seabed depths ranging from 1 m to 10 m, and moderately turbid water conditions, corresponding to a chlorophyll-a concentration C h l = 1 mg m 3 , a suspended matter concentration S M = 1 mg L 1 , and a color dissolved organic matter concentration C D O M = 0.1 m 1 .
Remotesensing 12 02072 g001
Figure 2. Direct and diffuse water column attenuation (respectively k 1 and k 2 vectors); (a) Clear water ( C h l = 0.03 mg·m 3 , S M = 0.01 mg·L 1 , C D O M = 0.01 m 1 ); (b) Moderately turbid water ( C h l = 1 mg·m 3 , S M = 1 mg·L 1 , C D O M = 0.1 m 1 ).
Figure 2. Direct and diffuse water column attenuation (respectively k 1 and k 2 vectors); (a) Clear water ( C h l = 0.03 mg·m 3 , S M = 0.01 mg·L 1 , C D O M = 0.01 m 1 ); (b) Moderately turbid water ( C h l = 1 mg·m 3 , S M = 1 mg·L 1 , C D O M = 0.1 m 1 ).
Remotesensing 12 02072 g002
Figure 3. Spectral angle mapper, mean normalized R M S E errors on estimated spectra and estimated abundances using unmixing methods for various cases: when no water column is present ( N M F , no water), when the water column is present but the adjacency effects induced by the seabed and the water column are ignored ( W U M ), and when the water column and the adjacency effects are accounted for ( W A D J U M ), for clear and moderately turbid water conditions. K 1 , K 2 and δ are known.
Figure 3. Spectral angle mapper, mean normalized R M S E errors on estimated spectra and estimated abundances using unmixing methods for various cases: when no water column is present ( N M F , no water), when the water column is present but the adjacency effects induced by the seabed and the water column are ignored ( W U M ), and when the water column and the adjacency effects are accounted for ( W A D J U M ), for clear and moderately turbid water conditions. K 1 , K 2 and δ are known.
Remotesensing 12 02072 g003
Figure 4. Normalized R M S E for endmember spectra ( N S R M S E ) and abundances ( N A R M S E ), for five separate constant seabed depth areas (plain lines), and for the entire image composed of all five areas (dashed lines), obtained using biased values of K 1 , K 2 and δ induced by a lack of knowledge for the water quality, in the case of moderately turbid water and z = [ 1 , 2 , 3 , 4 , 5 ] m.
Figure 4. Normalized R M S E for endmember spectra ( N S R M S E ) and abundances ( N A R M S E ), for five separate constant seabed depth areas (plain lines), and for the entire image composed of all five areas (dashed lines), obtained using biased values of K 1 , K 2 and δ induced by a lack of knowledge for the water quality, in the case of moderately turbid water and z = [ 1 , 2 , 3 , 4 , 5 ] m.
Remotesensing 12 02072 g004
Figure 5. (a) Colored composition of the airborne reflectance data of the Porquerolles Island, acquired using the Hyspex hyperspectral sensor (with sub-surface and sun glint correction), with the trace of the ROV transect (black), and local study zone (red rectangle); (b) Bathymetry provided by the L i t t o 3 D data; (c) Color image of the seabed given by the Ifremer ROV, showing clear and dark sand, brown P o s i d o n i a , and C a u l e r p a (green); (d) Dark substratum, C a u l e r p a t a x i f o l i a and brown colonized P o s i d o n i a ; (e) Green P o s i d o n i a on clear sand; (f) Spectral library collected in the coastal zone of Porquerolles Island, used for simulations (the endmembers used for the seabed generation are marked with a dot).
Figure 5. (a) Colored composition of the airborne reflectance data of the Porquerolles Island, acquired using the Hyspex hyperspectral sensor (with sub-surface and sun glint correction), with the trace of the ROV transect (black), and local study zone (red rectangle); (b) Bathymetry provided by the L i t t o 3 D data; (c) Color image of the seabed given by the Ifremer ROV, showing clear and dark sand, brown P o s i d o n i a , and C a u l e r p a (green); (d) Dark substratum, C a u l e r p a t a x i f o l i a and brown colonized P o s i d o n i a ; (e) Green P o s i d o n i a on clear sand; (f) Spectral library collected in the coastal zone of Porquerolles Island, used for simulations (the endmembers used for the seabed generation are marked with a dot).
Remotesensing 12 02072 g005
Figure 6. Various spectra representative of sand, P o s i d o n i a and substratum, collected in the studied area.
Figure 6. Various spectra representative of sand, P o s i d o n i a and substratum, collected in the studied area.
Remotesensing 12 02072 g006
Figure 7. Part 1; (a) Abundances provided by the ground truth ( A G T ); (b) LS coefficients ( C e s t ); (c) Abundances obtained using WADJUM ( A e s t ).
Figure 7. Part 1; (a) Abundances provided by the ground truth ( A G T ); (b) LS coefficients ( C e s t ); (c) Abundances obtained using WADJUM ( A e s t ).
Remotesensing 12 02072 g007
Figure 8. Part 2; (a) Abundances provided by the ground truth ( A G T ); (b) LS coefficients ( C e s t ); (c) Abundances obtained using WADJUM ( A e s t ).
Figure 8. Part 2; (a) Abundances provided by the ground truth ( A G T ); (b) LS coefficients ( C e s t ); (c) Abundances obtained using WADJUM ( A e s t ).
Remotesensing 12 02072 g008
Figure 9. (a) Coefficients ( C e s t ), estimated using LS inversion; (b) Airborne reflectance data of the Porquerolles Island, taken using the Hyspex hyperspectral sensor (sub-surface correction); (c) Abundances ( A e s t , estimated using W A D J U M ); (d) Bathymetry provided by the L i t t o 3 D database. The depth is between 1 . 8 m (white) and 6 . 1 m (black).
Figure 9. (a) Coefficients ( C e s t ), estimated using LS inversion; (b) Airborne reflectance data of the Porquerolles Island, taken using the Hyspex hyperspectral sensor (sub-surface correction); (c) Abundances ( A e s t , estimated using W A D J U M ); (d) Bathymetry provided by the L i t t o 3 D database. The depth is between 1 . 8 m (white) and 6 . 1 m (black).
Remotesensing 12 02072 g009
Table 1. Notations used in the paper.
Table 1. Notations used in the paper.
δ Environment vector: contribution of the target pixel
reflectance to diffuse upward radiance
δ Environment parameter, δ = m e a n ( δ )
L u Sub-surface upward radiance
L u d p Radiance due to photons from the water column which have
not interacted with the seabed
L d i r Radiance due to photons directly transmitted by the water
column after interaction with the seabed
L d i f Radiance due to photons scattered after interaction with
the seabed
E t o t Incident downward flux that reaches the seabed
E 0 Incident downward flux that reaches sub-surface level
T d i r Direct upward transmittance of water column
T d i f Diffuse upward transmittance of water column
ρ t Reflectance of the target pixel
ρ e n v Equivalent reflectance of the target pixel environment
that contributes to the upward radiance, including the target
ρ a d j Equivalent reflectance of the pixels adjacent to the target
that contributes to the upward radiance, not including the
target
LNumber of wavelengths
INumber of pixels in the whole image
JNumber of endmembers
K 1 Direct attenuation matrix of upward light for the whole
image, with size L × I
K 2 Diffuse attenuation matrix of upward light for the whole
image, with size L × I
k 1 Direct attenuation vector for one pixel
k 2 Diffuse attenuation vector for one pixel
R Matrix of observed sub-surface reflectance, with size L × I
R ˜ Matrix R minus the water reflectance
S matrix containing endmember spectra, with size L × J
A matrix containing abundance coefficients, with size J × I
zseabed depth(bathymetry), in m
Table 2. Maximum relative spread of the total downward radiance E t o t at the seabed level, due to the variability of the seabed reflectance and coupling illumination.
Table 2. Maximum relative spread of the total downward radiance E t o t at the seabed level, due to the variability of the seabed reflectance and coupling illumination.
Water Column Condition c 1 c 2 c 3
m a x r ( E t o t ) 1 . 3 % 1 % 0 . 4 %
Table 3. Performance of W U M and W A D J U M when a random error with an uniform law U [ 0 . 5 , 0 . 5 ]  m is added to the seabed depth, for a true seabed depth of 5 . 5  m and of 9 . 5  m, in the case of moderately turbid water.
Table 3. Performance of W U M and W A D J U M when a random error with an uniform law U [ 0 . 5 , 0 . 5 ]  m is added to the seabed depth, for a true seabed depth of 5 . 5  m and of 9 . 5  m, in the case of moderately turbid water.
True Depth NK 1 RMSE NK 2 RMSE NS RMSE WUM WADJUM NA RMSE WUM WADJUM
5 . 5  m 6 % 5 % 22 % ± 7 % 4 % ± 1 % 30 % ± 11 % 14 % ± 2 %
9 . 5  m 48 % 18 % 41 % ± 3 % 22 % ± 7 % 60 % ± 7 % 29 % ± 12 %
Table 4. Normalized R M S E between the true K 1 , t r u e , K 2 , t r u e and the biased K 1 , K 2 , due to an erroneous knowledge for water quality, for  z = [ 1 , 2 , 3 , 4 , 5 ] m , and moderately turbid water.
Table 4. Normalized R M S E between the true K 1 , t r u e , K 2 , t r u e and the biased K 1 , K 2 , due to an erroneous knowledge for water quality, for  z = [ 1 , 2 , 3 , 4 , 5 ] m , and moderately turbid water.
z ( m ) 12345
N K 1 R M S E 5 % 14 % 32 % 68 % 138 %
N K 2 R M S E 4 % 4 % 3 % 3 % 3 %

Share and Cite

MDPI and ACS Style

Guillaume, M.; Minghelli, A.; Deville, Y.; Chami, M.; Juste, L.; Lenot, X.; Lafrance, B.; Jay, S.; Briottet, X.; Serfaty, V. Mapping Benthic Habitats by Extending Non-Negative Matrix Factorization to Address the Water Column and Seabed Adjacency Effects. Remote Sens. 2020, 12, 2072. https://doi.org/10.3390/rs12132072

AMA Style

Guillaume M, Minghelli A, Deville Y, Chami M, Juste L, Lenot X, Lafrance B, Jay S, Briottet X, Serfaty V. Mapping Benthic Habitats by Extending Non-Negative Matrix Factorization to Address the Water Column and Seabed Adjacency Effects. Remote Sensing. 2020; 12(13):2072. https://doi.org/10.3390/rs12132072

Chicago/Turabian Style

Guillaume, Mireille, Audrey Minghelli, Yannick Deville, Malik Chami, Louis Juste, Xavier Lenot, Bruno Lafrance, Sylvain Jay, Xavier Briottet, and Veronique Serfaty. 2020. "Mapping Benthic Habitats by Extending Non-Negative Matrix Factorization to Address the Water Column and Seabed Adjacency Effects" Remote Sensing 12, no. 13: 2072. https://doi.org/10.3390/rs12132072

APA Style

Guillaume, M., Minghelli, A., Deville, Y., Chami, M., Juste, L., Lenot, X., Lafrance, B., Jay, S., Briottet, X., & Serfaty, V. (2020). Mapping Benthic Habitats by Extending Non-Negative Matrix Factorization to Address the Water Column and Seabed Adjacency Effects. Remote Sensing, 12(13), 2072. https://doi.org/10.3390/rs12132072

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop