Spectral unmixing of Multispectral Lidar signals

In this paper, we present a Bayesian approach for spectral unmixing of multispectral Lidar (MSL) data associated with surface reflection from targeted surfaces composed of several known materials. The problem addressed is the estimation of the positions and area distribution of each material. In the Bayesian framework, appropriate prior distributions are assigned to the unknown model parameters and a Markov chain Monte Carlo method is used to sample the resulting posterior distribution. The performance of the proposed algorithm is evaluated using synthetic MSL signals, for which single and multi-layered models are derived. To evaluate the expected estimation performance associated with MSL signal analysis, a Cramer-Rao lower bound associated with model considered is also derived, and compared with the experimental data. Both the theoretical lower bound and the experimental analysis will be of primary assistance in future instrument design.


I. INTRODUCTION
L ASER altimetry (or Lidar) is an acknowledged tool for extracting spatial structures from three-dimensional (3D) scenes, including forest canopies [1], [2]. Using time-of-flight to create a distance profile, signal analysis can recover tree and canopy heights, leaf area indices (LAIs) and ground slope by analyzing the reflected photons from a target. Conversely, passive multispectral (MSI) (dozen of wavelengths) and hyperspectral images (HSI) (hundreds of wavelengths) are widely used to extract spectral information about the scene which, for forest monitoring, can also provide useful parameters about the canopy composition and health (tree species, leaf chlorophyll content, water content, stress, among others) [3], [4]. The most natural evolution to extract spatial and spectral information from sensed scenes is to couple Lidar data and multi/hyperspectral images [5], [6]. Although the fusion of Lidar data and HSIs can improve scene characterization, there are problems of data synchronization in space (alignment, resolution) and time (dynamic scene, change of observation conditions, etc). For these reasons, multispectral Lidar (MSL) has recently received attention from the remote sensing community for its ability to extract both structural and spectral information from 3D scenes [7], [8].
The key advantage of MSL is the ability to provide information on the vertical distribution of spectra, used to infer physiological processes directly linked to actual carbon sequestration as well as carbon stocks [9]. For example, a key capability is to directly measure and classify ground-based shrub infestation, which is difficult with conventional HSI as the lower spectral response is obscured by the 'first' signals returned from the top of the canopy. Preliminary trials with existing commercial LiDAR systems with three operational wavelengths have already taken place [10], but as we shall demonstrate here, such systems cannot yet provide all the necessary spectral and structural discrimination to directly measure these processes.
Another motivation for MSL is that HSI, even when fully synchronised, can only integrate the spectral response along the path of each optical ray, not measure the spectral response as a function of distance, e.g., depth into a forest canopy. Multiple scattering effects cannot be neglected in some scenes with relief or containing multi-layered objects (such as trees), which complicates surface detection and quantification. Most existing spectral unmixing (SU) techniques rely on a linear mixture assumption [11]- [15] to identify the components (so-called endmembers) of an image and their proportions (abundances). However, it has been shown that the classical linear mixing model (LMM) can be inaccurate when multiple scattering effects occur [16], [17]. Although polynomial models [18], [19] (including bilinear models [20]- [23]) can be adapted for "long range" multiple scattering (in contrast to intimate mixtures models [24]), motivating the additional parameter constraints and designing accurate nonlinear unmixing methods are still challenging problems [25]. Since MSL data present an additional dimension when compared to HSIs, we can expect a reduction of the SU problem complexity and thus an improvement in the target characterization performance.
In [26], the authors proposed a Bayesian algorithm to estimate the positions and amplitudes in Lidar signals associated with a multi-layer target. This method has been extended in [9], [27] to MSL by first estimating the positions of the peaks (i.e., the layers), which were assumed to be the same in all spectral bands, then estimating the amplitudes of the peaks for each wavelength, and finally relating the estimated amplitude peaks to areas associated with each material that make up the target based on a linear mixing model. The method has been applied to real MSL signals (four wavelengths) to analyze a conifer structure and has been shown to aid the recovery of bark and needle areas of the tree assuming it is modeled as a set of irregularly spaced layers.
In this paper, we propose to investigate a SU problem applied to single-and multi-layered targets to estimate the positions and proportions of each (known) material. This problem extends the This work is licensed under a Creative Commons Attribution 3.0 License. For more information, see http://creativecommons.org/licenses/by/3.0/ supervised SU problem of HSIs by also estimating the target position. Estimating the spectral variation of material signatures (unsupervised SU problem) is also of interest. For example, the estimation of physiological content, such as Chlorophyll concentration, in forest canopies is of key concern, and this results in variable spectra for leaves and needles. Extracting these spectra from MSL signals is a more challenging problem (probably more difficult than in HSIs due to the statistical properties of the noise and observation model) which is out of scope of this paper, but will require particular attention in future work.
In contrast with the classical additive Gaussian noise assumption used for HSIs, a Poisson noise model is more appropriate for MSL signals. Indeed, Lidar and thus MSL systems usually record, for each pixel/region of the scene, a histogram of time delays between emitted laser pulses and the detected photon arrivals. Within each histogram bin, the number of detected photons follows a discrete distribution which can be approximated by a Poisson distribution due to the particle nature of light. Using a Bayesian approach, appropriate prior distributions are chosen for the unknown parameters of the model considered here, i.e., the material areas, the surface positions and the background parameters. The joint posterior distribution of these parameters is then derived. Since the classical Bayesian estimators cannot be easily computed from this joint posterior (mainly due to the model complexity), a Markov chain Monte Carlo (MCMC) method is used to generate samples according to the posterior of interest.
More precisely, following the principles of the Gibbs sampler, samples are generated according to the conditional distributions of the posterior. Due to the possibly high correlations between the material proportions/areas, we propose to use a Hamiltonian Monte Carlo (HMC) [28] method to sample according to some of the conditional distributions. HMCs are powerful simulation strategies based on Hamiltonian dynamics which can improve the convergence and mixing properties of classical MCMC methods (such as the Gibbs sampler and the Metropolis-Hastings algorithm) [29], [30]. These methods have received growing interest in many applications, especially when the number of parameters to be estimated is large [31], [32]. Classical HMC can only be used for unconstrained variables. However, new HMC methods have been recently proposed to handle constrained variables [29,Chap. 5], [33], [34] which allow HMCs to sample according to the posterior of the Bayesian model proposed for SU. Finally, as in any MCMC method, the generated samples are used to compute minimum mean square error (MMSE) estimators as well as measures of uncertainties such as confidence intervals.
Predicting the parameter estimation performance is of prime interest for designing an estimation procedure but also assists with the instrument design. Indeed, since MSL is a recent modality, it is important to identify the key parameters that have an influence on the material estimation performance (such as the necessary signal to background levels, the number of bands to be considered and the associated wavelengths). To guide future instrument design, we consider a Cramer-Rao lower bound (CRLB) associated with the observation model and show that it can be used to approximate the estimation errors of the proposed Bayesian algorithm.
The remainder of the paper is organized as follows. Section II introduces the observation model associated with MSL returns for a single-layered object to be analyzed. Section III presents the hierarchical Bayesian model associated with the proposed model and its posterior distribution. The hybrid MCMC method used to sample from the posterior of interest is detailed in Section IV. Section V investigates the expected SU performance using Cramer-Rao lower bounds. Experimental results are shown and discussed in Section VI. Conclusions are reported in Section VII.

II. PROBLEM FORMULATION
This section introduces the observation statistical model associated with MSL returns for a single-layered object which will be used in Sections III and IV to solve the SU problem of MSL data. Precisely, we consider a set of observed Lidar waveforms where is the number of temporal (corresponding to range) bins. To be precise, is the photon count within the th bin of the th spectral band considered. Let be the position of a complex object surface or surfaces at a given range from the sensor, composed of materials whose spectral signatures (observed at wavelengths) are denoted as . According to [26], each photon count is assumed to be drawn from the following Poisson distribution (1) where is the photon impulse response whose shape can differ between wavelength channels, denotes the area of the th material composing the object (relative to a reference area) and is a background and dark photon level, which is constant in all bins at a given wavelength. Without loss of generality the photon impulse responses are assumed to be the same and modeled by the following piece-wise exponential is a positive hyperparameter vector. In this study, we assume that the shape parameters of the model (which is appropriate for our own lidar sensors) are fixed and known from the instrumental response. As seen later in Fig. 3 the shape is slightly asymmetrical when compared to the more common Gaussian model [35] which is used for higher power, longer pulse duration lidar systems. The consideration of a piece-wise exponential approximation is motivated by the actual shape of our instrumental impulse response but it is worth mentioning that any positive analytical approximation, even different approximations for the different wavelengths, could be used instead without modifying the proposed Bayesian estimation procedure. The consideration of a single impulse response model is done only for ease of reading. While a Gaussian approximation can lead to a slight bias in depth estimation, we shall use a Gaussian approximation in Section V to compute the Cramer-Rao bounds since the Gaussian approximation is twice differentiable (in contrast to the piece-wise exponential approximation), which simplifies the derivation of the bounds.
Due to physical considerations, the relative areas are assumed to satisfy the following positivity constraints (2) Similarly, the average background level in each channel satisfies . Spectral unmixing of the MSL data consists of estimating the position (i.e., ) of the target, the area of each component of the target, as well as the background levels from the observed data in . The next section studies a Bayesian model to estimate the unknown parameters in (1) while ensuring the positivity constraints mentioned above.

III. BAYESIAN MODEL
The unknown parameter vector associated with (1) contains the relative areas , the position of the target and background levels (satisfying positivity constraints). This section summarizes the likelihood and the parameter priors associated with these parameters.

A. Likelihood
Assuming that the detected photon counts/noise realizations, conditioned on their mean in all bins and for all wavelengths, are independent, (1) leads to (3) where , and denotes the th row of .

B. Prior for the Target Position
Since we don't have prior information about the position of the target, the following uniform distribution (4) is assigned to . Note that the position is a real variable that is not restricted to the integer values in .

C. Prior for the Relative Areas
To reflect the lack of prior knowledge about the areas, the following truncated Gaussian prior (5) is assigned to each area , where denotes the Gaussian distribution restricted to , which hidden mean and variance 0 and , respectively. The hyperparameter is shared by all the parameters and is arbitrarily fixed to a large value to ensure a weakly informative prior. Assuming prior independence for the unknown areas yields (6) where means "proportional to" and is the indicator function defined on .

D. Prior for the Background Level
Similarly, assigning a truncated Gaussian prior to each background level and assuming prior independence between these parameters leads to (7) where the hyperparameter is shared by all the parameters and is arbitrarily fixed to a large value to ensure a weakly informative prior.
In this paper, we assume that prior knowledge about the relative areas and background levels is limited. We use weakly informative Gaussian priors restricted to the positive orthant to reflect this lack of knowledge and reduce the estimation bias while ensuring the positivity constraints are satisfied. However any proper weakly informative prior distributions could have been used, leading to similar estimation results. Similarly, if useful information about the unknown parameters is available, more informative (possibly physics-based) priors can be used instead of (4) and (6) (target position and composition) and (7) (background levels).

E. Joint Posterior Distribution
The resulting directed acyclic graph (DAG) associated with the proposed Bayesian model is depicted in Fig. 1. The joint posterior distribution of the unknown parameter vector can be computed using (8) where has been defined in (3) and . Unfortunately, it is difficult to obtain closed form expressions for standard Bayesian estimators associated with (8). In this paper, we propose to use efficient Markov Chain Monte Carlo (MCMC) methods to generate samples asymptotically distributed according to (8), employing a Gibbs sampler. The principle of the Gibbs sampler is to sample according to the conditional distributions of the posterior of interest [30,Chap. 10]. Due to the high correlation between the elements of (which will be discussed later in the paper), we use an Hamiltonian Monte Carlo (HMC) method which improves the mixing properties of the sampling procedure (compared to a classical Metropolis-within-Gibbs sampler using Gaussian random walks). This method is detailed in the next section.

IV. BAYESIAN INFERENCE USING A CONSTRAINED HAMILTONIAN MONTE CARLO METHOD
In this Section, we develop an HMC-within Gibbs sampler to generate samples according to (8) and estimate the unknown parameters involved in (1) in order to solve the SU problem of MSL data for a single-layered object. The proposed sampler consists of three steps to update sequentially and and is summarised in Algo. 1.

Algorithm 1: Gibbs Sampling Algorithm (single layer) 1 Fixed input parameters:
, number of burn-in iterations , total number of iterations 2 Initialization 3 Iterations 4 Sample from the pdf in (9) and CHMC 5 Sample from the pdf in (11) and a Gaussian random walk 6 Sample from the pdfs in (13) and Gaussian random walks 7 Set .

A. Sampling the Areas
The full conditional distribution of is given by (9) and is not a standard distribution which is easy to sample. Consequently, it is the norm to use an accept/reject procedure to update . A classical and simple approach would use a multivariate Gaussian random walk. However, in practice the elements of can be highly correlated, especially when the materials present similar spectral signatures (which will be the case for vegetation targets). Consequently, a Hamiltonian Monte Carlo method [29,Chap. 5] is preferred to improve the mixing properties of the sampler. The principle of HMCs is to introduce auxiliary, or momentum, variables, and perform a Metropolis-Hastings move in a higher dimensional parameter space. The proposal distribution is then built to take into account the shape of the target distribution (9). Due to the area constraints (2), a constrained HMC must be used. In this paper, we used a constrained HMC (CHMC) scheme similar to that described in [19] and thus consider the following potential energy function (10) to simulate Hamiltonian dynamics and compute the appropriate acceptance ratio (see [19], [29] for technical details). This choice of potential energy function ensures that where is a positive constant. Note that the proposed CHMC step can be applied for any differentiable relative area prior. Metropolis-Hastings or more complex HMC moves should be used instead when considering non-differentiable priors instead of (6).

B. Sampling the Target Position
It can be shown from (8) that (11) is not a standard distribution and an accept/reject procedure must be used to update the target position . We use a Gaussian random walk to update this parameter. More precisely, a doubly truncated Gaussian proposal is considered to ensure each candidate belongs to the admissible set and the variance of the proposal is adjusted during the burn-in period of the sampler to obtain an acceptance rate close to 0.45, as recommended in [36, p. 8].

C. Sampling the Background Levels
It can be shown from (8) that (12) where (13) Sampling from (13) is again not straightforward and Gaussian random walks are used to update the background levels, similar to the position . However the background levels are a posteriori independent and can be updated in parallel. Similar to the target position update, the variances of the parallel Gaussian random walk procedures are set during the burn-in period of the sampler to obtain an acceptance rate close to 0.45.
After generating samples using the procedures detailed above and removing iterations associated with the burn-in period of the sampler ( has been set from preliminary runs), the MMSE estimator of the unknown parameters can be approximated by computing the empirical averages of the remaining samples. The minimal length of the burn-in period can be determined using several convergence diagnostics [36] but the number of initial samples to be discarded generally varies between data sets and can be thus difficult to adjust in advance without overestimating it. In previous work on parallel acceleration of RJMCMC algorithms we have investigated the application of such diagnostics to our monochromatic LiDAR signals [37], but as absolute speed of processing is not a major concern here, the length of the burn-in period is assessed visually from the preliminary runs and fixed to ensure that the sampler has converged. The total number of iterations has then been set to obtain accurate approximations (computed over 4000 samples in Section VI) of the MMSE estimators.

V. PREDICTING UNMIXING PERFORMANCE
This section studies a Cramer-Rao lower bound associated with the observation model (1) which can be used to assess the performance of methods that aim to solve the SU problem of MSL data considered in this paper as well as to assist with the future instrument design. Precisely, Section V-A recalls the definition of the CRLB for the problem of interest and Section V-B discusses the impact of several key parameters on the expected parameter estimation performance of SU methods.

A. Cramer-Rao Lower Bound
Prediction of the unmixing performance is necessary to assist in the design of a lidar system for a specific application, identifying those parameters that have the most significant impact. In deriving a CRLB, we propose firstly to relax the impulse response model in (2) by considering the following Gaussian approximation (14) This simplifies the CRLB derivation and does not significantly bias the prediction as the piece-wise exponential impulse response can be accurately approximated by a Gaussian function. The CRLB associated with any unbiased estimator of the parameter vector involved in the mixing model (1) (having replaced by ) and constructed from is given by (15) where is the Fisher information matrix whose elements are 1 The th diagonal element of the CRLB matrix in (15) provides a lower bound for the variance of , given that is an unbiased estimator of . Of course, the Bayesian estimation procedure proposed in Sections III and IV does not provide a strictly unbiased estimator of and a Bayesian CRLB should have been considered instead of (15). However, due to the weakly informative priors chosen in Section III, when the actual vector is far enough from the boundaries of the admissible set (defined by the positivity constraints), the bias of the proposed estimation procedure can be neglected and the resulting estimator achieves the CRLB (as will be shown in Section VI). Moreover, the CRLB in (15) can also provide information about the estimation performance of a potential optimization method that could be proposed to estimate based on maximum likelihood estimation (MLE).

B. Performance Analysis
We first investigate the performance of spectral unmixing of MSL signals for a single layered, artificial target composed of materials, specifically needles, bark and soil. These materials have been chosen because of our interests in forest canopy monitoring using MSL signals [9]. The reflectance spectra of these materials, observed at equally spaced spectral bands ranging from 400 nm to 2500 nm are depicted in Fig. 2. It is interesting to note the strong similarity between the bark and needles spectral signatures, which can make their discrimination difficult and highlights the need of multiple wavelengths 1 The Fisher information matrix is derived in the Appendix  to quantify these materials. The maximum number of spectral bands considered in this paper has been set to , which is a realistic value for a short-term real measurement campaign (the current instrument uses only 4 wavelengths). The relative area of needles (resp. bark and soil) has been arbitrarily set to (resp. and ). These relative areas do not necessarily sum to one as we consider that part of the laser light can penetrate through the target (semi-transparent target or a significant gap fraction) and may not be further reflected (single hit assumption). The fixed model parameters have been fixed from experiments to (16) These parameters will be fixed in the remainder of the paper unless otherwise specified. The impulse response parameters have been set to , , , , , and for the piece-wise exponential approximation in (2) and for the Gaussian approximation in (14). These parameters have been obtained by fitting the experimental impulse response measured in [26]. Fig. 3 shows that the Gaussian approximation provides a good estimate of the experimental impulse response, although the piece-wise exponential better fits the experimental curve. Fig. 4 shows an example of MSL data generated using the parameters in (16) and the impulse response approximation in (2). Table I shows the predicted estimation performance for the unknown parameters of interest (i.e., and ) assuming the Gaussian approximation of the photon impulse response (14). The relative errors provided in the bottom row of Table I are computed by dividing the square root of the CRLBs by the actual values of the parameters. The relative estimation errors of the background levels (not presented here) are lower than 1%.
These results show that the estimation errors associated with the target position are usually much lower than those associ-   ated with the material areas. Moreover, although the amplitude of the peak for each wavelength is much larger than the background level (see Fig. 4), the CRLB predicts possibly large estimation errors, especially for the needle and bark areas. This can be partially explained by the fact that the soil reflectance spectrum presents an average energy which is higher than those of needles and bark (and thus ) but also and mainly because of the high correlation between the bark and needle spectra, which complicates their quantification. Fig. 5 to Fig. 7 show the predicted estimation errors for different values of the key parameters (which is related to the number of photons emitted by the laser sources), the number of spectral bands (equally spaced between 400 nm and 2500 nm) and the average background level (assumed to be the same in all spectral bands and bins). Fig. 5 shows that the estimation performance generally increases with the number of spectral bands. This result is well known for SU of multispectral and hyperspectral images and is here demonstrated for MSL signals. However, the performance improvement can vary, depending on the additional spectral bands considered. For instance, additional bands containing similar spectral information may not significantly improve the unmixing results (e.g., several wavelengths between 2000 nm and 2500 nm for the three materials in Fig. 2). Fig. 6 shows that the parameter has a significant impact on the SU performance. This parameter increases with the amplitude and the number of laser source pulses used to acquire the data and decreases with the average distance between the source and the target (as the probability of recording a reflected photon decreases). Increasing leads to better area estimation but can require a longer target exposure (which can be problematic for airborne sensors for instance). Finally, Fig. 7 shows that for relatively small values of background levels (compared to ), the SU performance is not significantly degraded when the background increases. This background level (which depends primarily on the background (e.g., solar) radiation as well as the instrument design) is expected to be quite small in practice.

VI. EXPERIMENTS
In this section, we first apply (Section VI-A) the proposed SU method to synthetic MSL data associated with a single-layered target and compare its estimation performance to that predicted by the CLRB presented in Section V and to that of a state-of-the art method [9]. In Section VI-B, we investigate the case of a multi-layered target whose layer positions are assumed to be known. Although the current algorithm does not estimate the positions of multiple surfaces, Section VI-B aims to illustrate how the proposed method can be extended to analyze more complex objects. The main steps of the methodology (proposed to analyze single and multi-layer target) are summarised in Fig. 8.

A. Single Layer Target
In this section, we evaluate the performance of the proposed spectral unmixing algorithm on synthetic MSL signals generated using the parameters used in Section V-B ((16)) and impulse response approximation (2).
The number of spectral bands has been chosen as , equally spaced between 400 nm and 2500 nm. For each scenario, the number of iterations has been set to (including burn-in samples). The  hyperparameters have been fixed to . The estimation performance of the proposed algorithm is evaluated by comparing the CRLB defined in Section V to the mean square error (MSE) defined as (17) where and are the actual and estimated th element of the unknown parameter vector and the expectation is approximated using 900 Monte Carlo runs. The performance of the proposed algorithm is compared to the CRLB (including the Gaussian approximation (14)) and to the performance of the algorithm developed in [9], denoted as "State-of-the-art" and which consists of estimating sequentially the position, amplitudes of the peaks for each wavelength and finally the relative areas. Table II shows that the MSEs obtained for by the two algorithms are close to the CLRBs, although the proposed algorithm generally outperforms the method proposed in [9] due to the joint estimation of the target position and the material areas. Note that for , the variance of the proposed estimator is slightly lower than the CRLB. This can be mainly explained by the fact that the estimator is no longer unbiased in that case. Table III compares the average relative estimation errors obtained by the two algorithms and computed by dividing the square root of the MSEs by the actual values of the parameters. These results show that increasing the number of spectral bands from to almost divides the area estimation errors by three (e.g., from to for the bark area), which highlights the benefits of increasing the number of wavelengths.

B. Extension to a Multi-Layer Target
In this section we extend the model (1) by considering a target composed of layers, leading to (18) where is the photon impulse response of the th layer located at and is the relative area vector associated with the th layer. In this scenario, we assume that the number and positions of the layers are known. Consequently, the unmixing problem reduces to estimating the relative areas of the known components for the layers (and the background levels). The joint spectral and spatial unmixing problem is more challenging and will be discussed in Section VI-C. The Bayesian model presented in Section III has been extended by considering the same prior (6) for the area vectors of each layer. The sampling procedure studied in Section IV has been modified in order to update sequentially the area vectors and the target position update step has been removed since the surface positions are assumed to be fixed in this paragraph (see Algo. 2).

1) Target Description:
We evaluate the unmixing performance of the extended Bayesian algorithm using an artificial target composed of layers. More precisely, the multilayer target is assumed to be far enough from the source to ensure that the laser rays hitting the target have the same incident direction. Thus the photons reflected onto the th layers lead to the same impulse response (i.e., same ). The first two layers are located at and and are composed of needles, bar and soil. It is important to note that the 3 materials composing the first two layers are assumed to randomly distributed, without overlap. The third layer modeling the reference spectralon is located at . The associated relative areas are presented in Table IV and the location of the layers is depicted in Fig. 9. Note that in this particular scenario, the relative areas correspond to areas visible by the source (no occlusion) and satisfy . However, this constraint is not ensured by the proposed estimation procedure. MSL data associated with this scene have been generated according to (18) with and and are depicted in Fig. 10 2) Estimation Performance: The proposed Bayesian algorithm extended to account for multiple layers has been applied to the MSL signals with (including burn-in samples). The layer positions have been fixed to their actual values. Table IV shows the estimated relative areas for  each layer and . This table shows that the proposed algorithm provides accurate area estimates for each layer and the more wavelengths the better the estimation, as already observed with the single-layered target in Section VI-A.

C. Towards Real Data Analysis
Due to the current limitations of the instrument, the proposed algorithm has only been applied to synthetic MSL signals whose underlying model has been shown to be in good agreement with our previous real MSL measurements using fewer wavelengths [9], [26]. Using this real data, we have managed to investigate the precision of area estimation in a multi-layer model; ground truth was available for the endmembers, but unfortunately we did not have structural truth. As stated above, real MSL signal acquisition campaigns are underway with 3 wavelengths [10] but we need to extend these and have simultaneous ground truth measurements to better evaluate the accuracy of the model and the estimation algorithm presented in this paper. Moreover, the number and positions of the layers which the multi-layer target is composed of were assumed to be known in Section VI-B. Estimating these parameters, especially the number of layers, for MSL data is a challenging problem that can be addressed in the Bayesian framework using reversible jump MCMC methods, in a similar fashion to [26]. This issue is however out of scope of this work and will also be addressed in a future paper.

VII. CONCLUSION
In this paper, we have proposed and developed a Bayesian model and an MCMC method for spectral unmixing of multispectral Lidar data. First, we evaluated our method on a single layer simulated target composed of known materials so that the MSL returns consisted of the sum of the individual contributions of the different components. The algorithm estimated the target position, the relative areas of the materials and the noise statistical properties.
The Cramer-Rao lower bound associated with the observation model was derived to identify the key parameters involved in the expected performance of the SU algorithm. It was shown that this bound can be used to help design future multispectral/hyperspectral Lidar systems (e.g., number of spectral bands and associated wavelengths, laser power ...) for specific applications. It was also shown that the performance of the proposed SU strategy is close to the results provided by the Cramer-Rao lower bound, although the bound considered in the paper only holds for deterministic parameters. The proposed Bayesian model was then extended to handle multi-layer targets (assuming the layer positions are known) and the simulation results provided interesting results for solving the more challenging joint spectral and spatial unmixing problem.
The model use in this paper does not take into account possible multiple scattering effects (which are likely to occur when analyzing multi-layer scenes). Such effects have already been observed in HSIs and Lidar signals over canopies. Developing non-linear models for MSL signal analysis is a challenging problem to be addressed in future studies.

APPENDIX FISHER INFORMATION MATRIX
The likelihood of the observation matrix can be expressed as (19) where . The corresponding log-likelihood is given by (20) The partial derivatives of with respect to (w.r.t.) the unknown model parameters are Straightforward computations lead to if else using the fact fact