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

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 (NMF)-based unmixing algorithms, referred to as WADJUM (Water ADJacency UnMixing) and WUM (Water UnMixing, no adjacency effects) are developed, assuming as known the water column bio-optical properties. Simulations show that WADJUM algorithm achieves performance close to that of the NMF-based unmixing of the seabed without any water column, up to 10 m depth. WUM 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 WADJUM approach. WADJUM 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 (LS) estimation of the seabed mixing coefficients. The results show the relevance of the WADJUM algorithm for the local analysis of the benthic habitats.


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 OSOAA 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 (Chl), suspended non-algal particles (SM), and colored dissolved organic matter (CDOM) [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 MF) 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 MF [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 LIGU. A variant of LIGU is also used in the ZyCI AT toolbox by [35]. Seabed mapping is also developed in [36]. In [37], Maritorena's radiative transfer model [13] has been used to combine iteratively ML estimation and N MF 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 (ADMM)-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 MF-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 MF-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 MF-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 MF-based unmixing method, hereafter referred to as WADJU M (Water ADJacency UnMixing), is proposed. A simpler variant, referred to as WU 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 WADJU 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 LS 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 OSOAA 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.

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. 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 × Ĩ 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 z seabed depth(bathymetry), in m

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): (1) L dp u 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 dir and L di f are both related to the downward flux, to the seabed reflectance and to the water column transmittance: The symbol denotes the elementwise product and each bold term corresponds to a spectral vector. ρ t and ρ env respectively refer to the target (pixel P) reflectance and the environment (P and adjacent pixels) reflectance.
(a) (b) 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 Chl = 1 mg m −3 , a suspended matter concentration SM = 1 mg L −1 , and a color dissolved organic matter concentration CDOM = 0.1 m −1 .

Downward Flux E tot
The total incident downward flux that reaches the pixel P, E tot (P), is the sum of four components, corresponding respectively to the direct, diffuse, coupling and reflected light (see Figure 1a): The coupling term E coup 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 tot due to the multiple reflections between the seabed and the water column, described by the coupling term E coup . 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 maxr(R) = 61% for the seabed. The value of E tot has been generated using OSOAA 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 [Chl = 0.03 mg·m −3 , SM = 0.01 mg·L −1 , CDOM = 0.01 m −1 ], and the moderately turbid is defined as [Chl = 1 mg·m −3 , SM = 1 mg·L −1 , CDOM = 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 tot over the 2400 pixels have been calculated, and the resulting variability index, maxr(E tot ), is reported in Table 2. Table 2. Maximum relative spread of the total downward radiance E tot at the seabed level, due to the variability of the seabed reflectance and coupling illumination.

Water Column Condition
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 tot with the seabed composition can be assimilated to (weak) noise. Then multiple reflections that cause the variations of E tot can be neglected in the expression of the radiative transfer model of Equation (2). In the following of the paper E tot 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.

Environment Function and Environment Reflectance Analysis
In [45], the authors make use of the environment function G env , 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 env 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 ∑ L λ=1 δ ↑ (λ), 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: where ρ t is the target reflectance and ρ adj is the reflectance of the pixels adjacent to the target.
To build the underwater mixing model, it is necessary to decompose ρ t and ρ env as mixtures of endmember spectra.

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.
(ii) As mentioned previously in Section 2.2, OSOAA 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 tot , 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 di 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 : 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 ρ adj = 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: The 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 .

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 dp u (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: For each pixel i of the whole data, the following notations are introduced: By using the mixture decomposition of ρ env and ρ t , provided in Equation (6), the expression in (7) becomes:r For the entire scene containing I pixels and L wavelengths, the proposed sub-surface mixing model, denoted as SM wadj in the following, is then: The L × I matricesR, 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 SM wnoadj , is obtained:

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 LS inversion method that would not require the accurate knowledge of the seabed (for example SW I M code [48]), or alternatively by inversion using a fixed spectral library [36].
Then E 0− , E tot , L dp u , T dir , T di f , δ can be numerically calculated for each pixel, using the OSOAA 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 RQE, defined in Equation (12), is minimized.
The RQE is the Frobenius norm of the difference between the observed sub-surface reflectance modified according to (8) (see also Table 1),R obs , 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 MF are developed hereafter, respectively so-called WADJU M, corresponding to the mixing model SM wadj expressed in Equation (10), and WU M, corresponding to the mixing model SM wnoadj of Equation (11). The two variants are both based on the Minidisco 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 (STU) for the abundances, or the minimum dispersion constraint for the spectra [50]. The estimation of the matrices A, S is obtained following Equation (13): Λ denotes regularization terms that can be added to the RQE criterion. In our experiments, a soft constraint to encourage the sum-to-unity for abundances is used, The vector I 1J is composed of the values 'one'; the regularization coefficient λ STU is empirically set to 0.5 after various trials. For the WADJU M algorithm, the gradients of RQE, relative to A and S, are detailed in Equation (14), where t denotes the transpose. (14) In the WADJU M unmixing algorithm, updates of A, S, andR 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 ofR).
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 WU M, is obtained from the same algorithm, by setting P = Id I×I .

Analysis of Performance Based on Simulations
Simulations are carried out to evaluate the performance of the proposed WU M and WADJU M algorithms.

Evaluation Metric
The results are evaluated using the mean spectral angle mapper (SAM) and the mean normalized RMSE, defined as follows in Equation (15): Here ω and i respectively refer to the realizations (different noise generation and random initialization) and to the pixels, s i andŝ i refer to the true and estimated spectra, respectively; M andM 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 RMSE.
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 WU M and WADJU M are compared to those of the simple unmixing of the seabed with no water column, using the algorithm Minidisco [50], hereafter referred to as 'N MF, 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 WADJU M and WU M.

Simulated Data
The influence of the water column optical properties (attenuation due to scattering and absorption by hydrosols) is accounted for using OSOAA model simulations, using the selected water turbidity and depth as inputs of OSOAA. 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 dp u , E tot , E 0− , T dir , T di f , δ ↑ were calculated using the OSOAA 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 reflectanceR 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 tot at the seabed level and the direct transmittance T dir both systematically decrease with depth, whereas the diffuse transmittance T di 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 di f E tot π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.

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 WADJU M algorithm (SM wadj model), and the WU M algorithm (SM wnoadj model). The unmixing method is also applied in the absence of a water column, using the Minidisco algorithm [50], denoted as 'N MF, no water' in the figures. The results of Minidisco 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 WADJU M and WU 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 FCLS (Fully Constrained Least-Square) algorithm [52] applied to the seabed. The algorithms were developed in the MATLAB language. Minidisco (N MF, 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 WADJU M and WU 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 (Chl = 0.03 mg·m −3 , SM = 0.01 mg·L −1 , CDOM = 0.01 m −1 ), and moderately turbid water (Chl = 1 mg·m −3 , SM = 1 mg·L −1 , CDOM = 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 (SAM), the normalized root mean square error (NS RMSE ) on estimated endmember spectra, and the N A RMSE 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. . Spectral angle mapper, mean normalized RMSE errors on estimated spectra and estimated abundances using unmixing methods for various cases: when no water column is present (N MF, no water), when the water column is present but the adjacency effects induced by the seabed and the water column are ignored (WUM), and when the water column and the adjacency effects are accounted for (WADJUM), for clear and moderately turbid water conditions. K 1 , K 2 and δ are known.
The results of the unmixing method using Minidisco, applied to the case where the water column is not present, show an angle error around 0.02 rad, a normalized RMSE 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 RMSE is around 6% for spectra, and the N A RMSE is around 12% for abundances. In clear water, the difference between WU M and WADJU 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 WU 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 WADJU M remains about the same as the one obtained in clear water context. Standard deviations are comparable for the WU M and WADJU M algorithms in the case of clear water, whereas the RMSE obtained using the WU M algorithm appears to be unstable for moderately turbid waters at the intermediate seabed depth of 5 m. Therefore, the proposed WADJU 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.

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 WU M and WADJU 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.

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 true = 5.5 m and z true = 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 OSOAA 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. Table 3. Performance of WUM and WADJUM 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.
For the deepest seabed case (z true = 9.5 m), the random error on the depth leads to a significant error on the attenuation matrices, NK1 RMSE = 47% and NK2 RMSE = 18%. This high error on the water attenuation results in degraded performance for both WU M (NS RMSE = 41%, N A RMSE = 60%) and WADJU M (NS RMSE = 22%, N A RMSE = 29%), although the estimation remains better for WADJU M.
In any case, the WADJU M algorithm still provides better results than the WU M algorithm, so the relevance of taking in consideration the adjacency effects for the unmixing is confirmed. In moderately turbid water, the WADJU 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.

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 WU M and WADJU 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 RMSE 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, NK1 RMSE , becomes significant when the depth increases, while NK2 RMSE 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, NS RSME , also increases with depth, although it is lower using WADJU M than WU 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, NS RSME , 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 RQE, 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 WU M and WADJU 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. WADJU M is more performing than WU 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 RMSE ) when using WADJU M in the case of biased water attenuation knowledge. The difference between WADJU M and WU 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.

Application to Real Data
An airborne data campaign has been conducted by the Hytech Imaging company [53] on 13 September 2017, using the HySpex 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 imageR (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 (Chl), suspended matter (SM) and colored dissolved organic matter (CDOM) 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 RGB 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 ROV 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 Chl = 0.03 mg·m −3 , SM = 0.5 mg·L −1 , CDOM = 0.05 m −1 . The bathymetry was provided by the LITTO3D database, available in [58], supported by the IGN (French National Institute for Geographical and Forestry Information) and the SHOM (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, Chl, SM, CDOM, and the bathymetry.

Initialization of the Unmixing
Initial matrices S init and A init 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 LS method (lsqcurvefit function of MATLAB), using a chosen library S 0 , and the model SM wnoadj of Equation (11), to obtain the coefficient matrix C est . 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 inv = S 0 * C est , is then recovered.
Step two : the initial matrix S init is obtained by applying the VCA algorithm [59] to X inv , and the initial abundance matrix A init is obtained using the FCLS algorithm [52] and S init .

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 est ) and coefficients (C est ) are shown, the estimated spectra being distorted due to the approximative knowledge of the water attenuation.

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 GT , the LS estimated coefficients C est , and the WADJU M estimated abundances A est are respectively shown in Figure 7a It can be seen that for both LS 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 WADJU M better fit the ground truth than the LS coefficients, which shows the interest of the unmixing method. The abundance errors obtained for the three unmixing methods are N A RMSE = 22% using VCA + FCLS performed on X inv , N A RMSE = 19% using WU M or WADJU M performed onR. As two classes are present, the error is the same for the sand and for the Posidonia. The error is not calculated for the coefficients C est , because they are not normalized.

Abundance Estimation Performed on Part 2
In the Part 2 of the validation area there is a Caulerpa 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 Caulerpa location. As the monitoring of the Caulerpa is an important issue for the PNPC, we focus on the estimation of such a benthic habitat. Posidonia 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 Caulerpa given by the LS coefficients C est is more dispersed than the distribution given by A est , obtained using WADJU M. The difference between the initialization using VCA + FCLS performed on X inv , and the output of WU M and WADJU M performed onR, is not tremendous for the estimation of Caulerpa abundance: N A RMSE = 15% using VCA + FCLS, N A RMSE = 13% using WU M, N A RMSE = 12% using WADJU M.
In conclusion, the validation shows a real improvement using an unmixing strategy instead of LS coefficient estimation, and a slight improvement using WU M or WADJU M performed on the sub-surface images compared to the use of VCA + FCLS performed on the reconstructed seabed with no water column. The difference between WU M and WADJU 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 VCA + FCLS and WU 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 WU M and WADJU M unmixing methods (see Section 3.4.2).

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 LS inversion were dark green Posidonia, 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 dp u 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 MATLAB, the image has been split into five contiguous sub-images to run the WADJU 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 Posidonia 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 WADJU 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 Posidonia seagrass is observed in the middle of the area identified as substratum/brown Posidonia, and photophilic algae are located on the edge of the area of settlement of the Posidonia 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 est (Figure 9a) and the abundances obtained in A est (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 Posidonia is erroneously detected on the seabed right corner of the coefficients in C est , and the substratum/brown Posidonia could be confused with sand in C est .

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 OSOAA 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 WADJU 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 WU 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 WU M and WADJU 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 WADJU M and WU 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 WADJU 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 ROV, in a depth range of 5.5 m-14 m. For these real data, the clear water condition induced no significant difference between WU M and WADJU 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 WU M and WADJU M were similar up to 10 m depth, for abundance estimation. Furthermore, the imperfect knowledge of the water quality also reduces the difference between WU M and WADJU M performances. On the other hand, the comparison of the estimated abundances A est with the coefficients C est 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 LS coefficient estimation using a fixed spectral library. In the validation area, the comparison with initialization matrices showed a small improvement using WADJU M, which could be due to imperfect knowledge of the water bio-optical properties, and also to the clear water conditions.

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 WU M algorithm can be used, while for moderately turbid waters and pixel sizes lower than 2 m, WADJU 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 (ADMM) 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 WADJU 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.

Conflicts of Interest:
The authors declare no conflict of interest.