Yingying Song

El-Hadi Djermoune

Jie Chen
email: dr.jie.chen@ieee.org

C Édric

David Brie

ONLINE DECONVOLUTION FOR INDUSTRIAL HYPERSPECTRAL IMAGING SYSTEMS *

Keywords: hyperspectral image, online deconvolution, LMS, ZA-LMS AMS subject classifications. 68U10

This paper proposes a hyperspectral image deconvolution algorithm for the online restoration of hyperspectral images as provided by wiskbroom and pushbroom scanning systems. We introduce a least-mean-squares (LMS)-based framework accounting for the convolution kernel non-causality and including non-quadratic (zero attracting and piece-wise constant) regularization terms. This results in the so-called sliding block regularized LMS (SBR-LMS) which maintains a linear complexity compatible with real-time processing in industrial applications. A model for the algorithm mean and mean-squares transient behavior is derived and the stability condition is studied. Experiments are conducted to assess the role of each hyper-parameter. A key feature of the proposed SBR-LMS is that it outperforms standard approaches in low SNR scenarios such as ultra-fast scanning.

Introduction. Hyperspectral imaging has received considerable attention in

the last decade as it combines the power of digital imaging and spectroscopy. Every pixel in a hyperspectral image provides local spectral information about a scene of interest across a large number of contiguous bands. This information can be used to characterize objects with great precision and details in a number of areas, including agricultural monitoring, industrial inspection, and defense. The core characteristics of hyperspectral images raise new data processing issues ranging from image restoration to pattern recognition [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF][START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF][START_REF] Ammanouil | Blind and fully constrained unmixing of hyperspectral images[END_REF][START_REF] Song | Regularization parameter estimation for non-negative hyperspectral image deconvolution[END_REF]. Several sensing techniques have been devised for hyperspectral imaging. They can be categorized into four main groups [START_REF] Willett | Sparsity and structure in hyperspectral imaging: sensing, reconstruction, and target detection[END_REF][START_REF] Li | Review of spectral imaging technology in biomedical engineering: achievements and challenges[END_REF]: whiskbroom (point scan), pushbroom (line scan), tunable filter (wavelength scan), and snapshot. Pushbroom systems are used in many areas such as food safety [START_REF] Chen | Transportable spectrophotometer system for on-line classification of poultry carcasses[END_REF][START_REF] Huang | Recent developments in hyperspectral imaging for assessment of food quality and safety[END_REF], georeferencing [START_REF] Cariou | Automatic georeferencing of airborne pushbroom scanner images with missing ancillary data using mutual information[END_REF] and material sorting [START_REF] Tatzer | Industrial application for inline material sorting using hyperspectral imaging in the NIR range[END_REF][START_REF]Pellenc Selective Technology[END_REF]. The main feature of pushbroom imaging systems is that the hyperspectral data cubes are acquired slice by slice, sequentially in time (which in fact represents the so-called along track spatial dimension). This paper is a first step towards the development of advanced online (real-time) hyperspectral image processing methods required in industrial processes that aim at controlling and sorting input materials right after each line scanning. For example, this is the case for paper and chipboard plants that use recycled materials in the manufacturing process. The aim of this paper is to address the fast online (sequential) deconvolution of hyperspectral images captured by pushbroom imaging systems. However, it is worth noting that the proposed sequential deconvolution algorithm can be easily extended to whiskbroom (point-by-point) systems.

Multichannel image restoration was carried out with Wiener methods in [START_REF] Hunt | Karhunen-Loeve multispectral image restoration, part I: Theory[END_REF][START_REF] Galatsanos | Digital restoration of multichannel images[END_REF].

Other strategies such as those in [START_REF] Galatsanos | Least squares restoration of multichannel images[END_REF][START_REF] Giovannelli | Positive deconvolution for superimposed extended source and point sources[END_REF][START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF][START_REF] Zhao | Deblurring and sparse unmixing for hyperspectral images[END_REF] were also introduced, but only in an offline setting. Let us mention the works of Jemec et al. [START_REF] Jemec | Push-broom hyperspectral image calibration and enhancement by 2D deconvolution with a variant response function estimate[END_REF][START_REF] Jemec | Deconvolution-based restoration of SWIR pushbroom imaging spectrometer images[END_REF] which specifically address the calibration and (offline) deconvolution of pushbroom hyperspectral images. The idea of restoring images in a sequential way can be traced back to the end of the 70's and has suscited extensive works up to the mid 90's. Early works on 2D-image restoration are due to Woods and Radewan, who developed a Kalman-filter (KF)-based approach allowing a pixel by pixel image restoration [START_REF] Woods | Kalman filtering in two dimensions[END_REF]. This approach was then modified in [START_REF] Woods | Kalman filtering in two dimensions: Further results[END_REF] to lower the computational cost, resulting in the so-called reduced update KF which was then extended to the restoration of color images [START_REF] Tekalp | Multichannel image modeling and kalman filtering for multispectral image restoration[END_REF]. An alternative KF-based approach was proposed by Biemond et al. [START_REF] Biemond | A fast Kalman filter for images degraded by both blur and noise[END_REF] to perform a lineby-line image restoration and which, due to Toeplitz matrix structures in the image and observation models, can be efficiently implemented in the Fourier domain. This approach was then extended to multichannel image restoration [START_REF] Galatsanos | Restoration of color images by multichannel Kalman filtering[END_REF] allowing a slice by slice restoration. The interest of KF-based image restoration started to decrease in the mid 90's with the growth of the computational capabilities. This led to the development of new iterative off-line image restoration approaches handling efficiently non-quadratic terms (including edge-preserving regularization) as in [START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. There is currently a renewed interest in extending these approaches to online estimation as [START_REF] Suzuki | Dual averaging and proximal gradient descent for online alternating direction multiplier method[END_REF][START_REF] Gao | Online dictionary learning for kernel LMS[END_REF][START_REF] Combettes | Stochastic forward-backward and primal-dual approximation algorithms with application to online image restoration[END_REF][START_REF] Chouzenoux | A stochastic majorize-minimize subspace algorithm for online penalized least squares estimation[END_REF] that comes from the need for processing tools that are able to handle very large data sets. However, extending such approaches to online hyperspectral image deconvolution can be tricky. For example, the online ADMM algorithm proposed in [START_REF] Suzuki | Dual averaging and proximal gradient descent for online alternating direction multiplier method[END_REF] addresses the problem of estimating fixed parameters and derives online algorithms that converge to the same solution as the batch ADMM. Basically they consist of (sub)-gradient-based algorithms with a step size that decreases to zero as the number of time samples increases. This may be irrelevant for online hyperspectral image deconvolution since, in that case, the key feature is the tracking capability of the algorithms.

The approach proposed in this paper has connections with [START_REF] Galatsanos | Restoration of color images by multichannel Kalman filtering[END_REF] since it addresses the slice-by-slice hyperspectral image restoration problem. It also uses the same approach to transform the non-causal blurring model into a (semi-)causal model. However, the KF equations need to update the state vector, the associated Kalman gain and covariance matrices. This complexity may not be compatible with real-time processing in industrial applications of hyperspectral imaging systems. The main motivation is to derive algorithms that allow to incorporate additional constraints enforcing some prior properties in the restored images while maintaining a linear complexity with respect to the number of unknown parameters. This work introduces an least mean-squares (LMS) framework for sequential deconvolution of hyperspectral images in which the subgradient is substituted by an instantaneous stochastic approximation and the step size is constant. Our algorithm operates in the spirit of LMS-based algorithms used, for instance, for adaptive system identification [START_REF] Widrow | Adaptive signal processing[END_REF][START_REF] Chen | Regularized least-mean-square algorithms[END_REF][START_REF] Chen | Nonnegative least-mean-square algorithm[END_REF] and superresolution restoration of sequences of images [START_REF] Elad | Superresolution restoration of an image sequence: adaptive filtering approach[END_REF][START_REF] Costa | Statistical analysis of the LMS algorithm applied to super-resolution image reconstruction[END_REF]. However, LMS-based sequential deconvolution of (hyperspectral) images has never been reported in the literature. Accounting for the specificities of the online image deconvolution (non-causal blurring, low signal-to-noise ratio), we propose a sliding-block regularized least mean-squares (SBR-LMS) algorithm that allows the restoration of images within a preset delay Q.

The paper is organized as follows. In Section 2, the image-blurring model is stated.

Causality of the convolution kernel and of the estimator are discussed. In Section 3, the hyperspectral image deconvolution problem is introduced. The SBR-LMS algo-rithm is proposed to address this problem. The criterion includes a spatial regularizer promoting the restoration of piecewise constant objects, a zero-attracting regularizer promoting the removal of the conveyor background, and a spectral regularizer promoting the smoothness of the spectral response. In Section 4, the transient behavior of the algorithm is analyzed in the mean and mean-squares sense. Its stability is also discussed. In Section 5, simulation results are presented to assess the accuracy of the convergence behavior model in both stationary and non-stationary scenarios. The role of each hyper-parameter is discussed. Numerical results are reported to assess the performance of the proposed approach and to compare it with state-of-the-art methods. Finally, the algorithm is applied to real hyperspectral images acquired under different experimental conditions.

Blurring and causality issues.

2.1. Scanning technologies and data structure. Recording a 3D hyperspectral data cube can be performed in several ways, depending on how the 3D signal is acquired, as a sequence of either 1D or 2D signals. Each configuration results in a specific imaging system. Here, we restrict our attention to whiskbroom and pushbroom scanners which are the most usual configurations in industrial applications.

Figure 1 depicts schematically how a hyperspectral image is captured by whiskbroom and pushbroom imagers. In both cases, the scene is scanned by moving the imager or its field-of-view (FOV) across the scene. For example, for material analysis, the samples to be imaged are carried by a conveyor moving at constant speed, while in remote sensing applications, the imaging system is carried on an airborne platform moves over the scene. In a whiskbroom imaging system, an 1D sensor (i.e., a spectrometer) is used to form the 3D data cube. Thus, the whole spatial scene is observed point-by-point yielding a stream of spectra which are stacked into a 3D data cube.

A mechanized angular movement using a scanning mirror sweeping allows to cover a spatial line from one edge to the other. Alternatively, pushbroom imaging systems make use of 2D sensors allowing to observe the scene line-by-line at each time instant.

The stream of spatial-spectral arrays is stacked to form the 3D-data cube. Following [START_REF] Schaepman | Imaging spectrometers, The SAGE handbook of remote sensing[END_REF], respective pros and cons of whiskbroom and pushbroom imaging systems are as follows. Because whiskbroom systems rely on 1D detector, the calibration is much simpler than with pushbroom systems. The disadvantages of this design include the presence of a mechanical scanning system, the shorter integration time than in pushbroom systems, and the image forming geometry which is dependent on the scanning speed, the scan mirror arrangement and the conveyor movement.

As compared to whiskbroom scan, pushbroom scanners allow for longer integration times for individual detector elements; the pixel dwell time is approximatively the inverse of the line frequency. In addition, there are distinct but fixed geometric relations between the pixels within a scan line. Since 2D sensors are used as focal planes in these systems, the uniform calibration of the detector response is critical. Note that in a combined analysis of signal-to-noise ratio (SNR), uniformity, and stability, pushbroom scanners might not necessarily outperform whiskbroom systems even though they have a longer integration time. For example, [START_REF]Pellenc Selective Technology[END_REF] developed a whiskbroom hyperspectral imaging system for ultra fast material sorting application allowing to reach up to 3 meters per second scanning velocity. However, pushbroom scanners are increasingly used in industrial applications.

A hyperspectral image is a 3D data cube whose dimensions are generally referred to as across-track, along-track and spectral dimensions. However, in the sequel, to avoid any confusion, they will be designated as spatial (across-track), time (along-track) and spectral dimensions. We will denote a hyperspectral image by smile effect, non-uniform illumination, and sensor response. These degradations may be corrected using standard calibration techniques [START_REF] Schaepman | Imaging spectrometers, The SAGE handbook of remote sensing[END_REF][START_REF] Esmonde-White | Minor distortions with major consequences: correcting distortions in imaging spectrographs[END_REF]. In this paper, we are primarily concerned with another type of distortion, namely, spatial distortions resulting in a loss of spatial resolution.

Y ∈ R N ×P ×K ,
In a whiskbroom imaging system, the spatial resolution is mainly controlled by the size of the beam footprint on the scanned object. Additional distortion can be caused by the conveyor motion but, in this paper, we shall assume that the integration time is small enough so that this distortion can be neglected. In general, the spatial sampling period ι s is chosen equal to the beam footprint size yielding a low resolution image. To increase the spatial resolution, it is necessary to reduce ι s , but this results in a spatial blurring of the hyperspectral image as sketched in Figure 3, where pixels are represented by the dots, and the beam footprint is represented by the circles.

The integration time T mainly acts on the noise level: a low value for T results in a low SNR. Consequently, assuming a constant acquisition velocity ι s /T , any resolution improvement results in an increase of both blurring and noise level. Assuming that the scanning angle is small, the beam footprint size can be considered as invariant with respect to the pixel position. We also assume that the footprint size is independent of the wavelength. Thus, the corresponding distortion can be modeled by linear invariant convolution. For pushbroom imaging systems, the quality of the acquired images can be significantly degraded by the blurring effect resulting from the miscalibration of the optical system [START_REF] Jemec | Push-broom hyperspectral image calibration and enhancement by 2D deconvolution with a variant response function estimate[END_REF][START_REF] Jemec | Deconvolution-based restoration of SWIR pushbroom imaging spectrometer images[END_REF][START_REF] Esmonde-White | Minor distortions with major consequences: correcting distortions in imaging spectrographs[END_REF]. According to [START_REF] Tao | A digital sensor simulator of the pushbroom offner hyperspectral imaging spectrometer[END_REF], there are many reasons that can cause spatial degradation of the pushbroom imaging system such as optical aberrations, slit size, and detector element size. In [START_REF] Jemec | Deconvolution-based restoration of SWIR pushbroom imaging spectrometer images[END_REF], the observed image corresponds to the 3D convolution of the original image with a spatial-spectral blurring filter varying with the spatial position. Here, the dependence of the spatial blurring with respect to the spatial position is neglected. We consider situations where the spectral response is smooth enough to choose a spectral sampling ι λ (controlled by the spectral binning) greater than the support of the spectral blurring and a spatial sampling ι s smaller than the support of the point spread function. This results in a degradation model involving only spatial blurring (possibly depending on the wavelength) of the hyperspectral image which is given by Equation (2.1) in the next section. The integration time is supposed to be small for fast scanning yielding a low SNR.

ι s K Convolution kernel (a) ι ′ s K ′ (b)
2.3. Causality. We shall now discuss issues related to the causality of the convolution kernel and associated estimates. Following [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF], hyperspectral image blurring can be seen as P simultaneous spatial convolutions. For each wavelength λ p , the blurred spatial image Y p ∈ R N ×K is given by the 2D convolution:

(2.1) where ⊤ denotes the transpose of a matrix. We use the same notation for the columns of X p . We assume a finite length blurring kernel of size L along the time dimension, centered around 0 which means that both past and future values of x p k contribute to the observation y p k . In order to make the blurring kernel causal

Y p = H ⋆p * X p + Z p where Y p ∈ R N
(2.3) ỹp col{ỹ p k } K k=1 .
where col{•} stacks its vector arguments on top of each other. The least squares criterion can be written as:

(2.4)

||ỹ p -F p x p || 2
where x p is built similarly to ỹp and F p is a Toeplitz-block-Toeplitz matrix of proper dimensions. To make the dependence of the criterion on x p k explicit, we introduce the following partitions:

x p = x p⊤ 1:k-1 , x p⊤ k , x p⊤ k+1:K ⊤ , (2.5)
F p = F p 1:k-1 , F p k , F p k+1:K (2.6)
where x p i:j col{x p k } j k=i and F p i:j is the the submatrix formed by the columns (i -1)N + 1 through jN of F p . This results in:

(2.7)

||ỹ p -F p x p || 2 = ||ỹ p -F p 1:k-1 x p 1:k-1 -F k+1:K x p k+1:K -F p k x p k || 2 .
It is now clear that optimally estimating x p k requires all the past estimations xp 1:k-1 and future estimations xp k+1:K , which precludes the derivation of a sequential estimator for x p k . To address these issues, we propose to produce the estimates xp k with a delay Q. More precisely, at time instant k, we shall estimate

x p k-Q+1 from Q observations (ỹ p k , ỹp k-1 . . . , ỹp k-Q+1 ), coarse posterior estimates (x p k , xp k-1 , . . . , xp k-Q+2
) refined as k increases, and anterior estimates (

x p k-Q , xp k-Q-1 , . . . , xp k-Q-L+2
) which are no longer updated. This is the key idea of the sliding-block LMS deconvolution algorithm developed in the next section.

3. Online image deconvolution.

3.1. Block Tikhonov. In this section, we introduce the Tikhonov-like approach proposed in [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] and a direct extension for online hyperspectral image restoration.

They will serve as reference methods to assess the performance of the proposed LMSbased method. This Tikhonov-like approach, originally developed in an offline setting, consists of seeking the minimum of a criterion composed of three terms: a data fitting term, a spatial regularizer, and a spectral regularizer:

(3.1) min x J (x) = 1 2 y -Hx 2 + η s 2 ∆ s x 2 + η λ 2 ∆ λ x 2 .
where x and y denote, respectively, the original and the observed vectorized hyperspectral images. Operator ∆ s corresponds to a Laplacian filter promoting the smoothness along both spatial and time dimensions, and ∆ λ corresponds to a firstorder derivative filter along the spectral dimension. Parameters η s and η λ are respectively the spatio-temporal and spectral regularization parameters. Matrix H is block-diagonal and corresponds to the convolution kernel. Because the criterion is quadratic, solution of problem (3.1) results in the following linear estimator:

(3.2) x = (H T H + η s ∆ T s ∆ s + η λ ∆ T λ ∆ λ ) -1 H T y.
Due to the block diagonal structure of H, the estimator (3.2) can be efficiently implemented in the frequency domain (see [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] for details).

The Block Tikhonov (BT) deconvolution approach addresses problem (3.1) in an online way by sequentially restoring spatio-spectral arrays X k ∈ R N ×P . Let x k and y k be the vectorized matrices X k ∈ R N ×P and Y k ∈ R N ×P , respectively:

(3.3) x k col {x p k } P p=1 , y k col {y p k } P p=1 .
where superscript p refers to the spectral band. This BT algorithm implements estimator (3.2) with slidding blocks of Q samples [x k-Q+1 , . . . , x k ], ∀k = Q, . . . , K and outputs the estimated vector xk-(Q-1)/2 as the final result3 . This algorithm works fast since only quadratic regularization terms are considered. However, when other constraints are needed, such as ℓ 1 -norm-based ones, no explicit solutions are available and time consuming iterative algorithms are required.

In the next section, we focus on extending the zero-attracting LMS (ZA-LMS) algorithm proposed in [START_REF] Chen | Regularized least-mean-square algorithms[END_REF] to online deconvolution, while accounting for some image specificities (non-causal blurring, presence of smooth regions separated by abrupt edges, and low SNR) resulting in the so-called Sliding-block regularized LMS (SBR-LMS).

-Q + 1, namely, x k , . . . , x k-Q+2 .
In what follows, to simplify notations, y k refers to the delayed observation ỹk andˆsymbol is omitted. This results in the following criterion:

J ( updated x k , . . . , x k-Q+1 | past estimates x k-Q , . . . , x k-Q-L+2 ) = Q q=1 E y k-q+1 - L ℓ=1 H ℓ x k-q-ℓ+2 2 + η z Q q=1 x k-q+1 1 + η s Q q=1 D s x k-q+1 1 + η λ Q q=1 D λ x k-q+1 2 (3.4)
where E{•} stands for the expectation operator,

• 1 = n |{•} n | denotes the ℓ 1 -
norm of its argument, and {•} n stands for the n-th entry of a vector. Matrix H ℓ blkdiag{H p ℓ } P p=1 is a block-diagonal matrix.

The regularizer ||D s x k-q+1 1 promotes the restoration of piecewise constant patterns along the spatial dimension. The first-order filtering operator D s is defined as:

(3.5)

D s I P ⊗ T N
where ⊗ stands for the Kronecker product, matrix I J denotes the J ×J identity matrix and T J is the Toeplitz matrix of size (J -1) × J with first column [1, 0, . . . , 0] and first row [1, -1, 0, . . . , 0]. The zero-attracting regularizer x p k-q+1 1 aims at removing the conveyor background. The choice of these regularization terms is thus mainly motivated by the targeted application, namely, the inspection of objects put on a conveyor belt. At a given wavelength, the response of the conveyor after background removal is close to zero while that of the objects is supposed to be piecewise constant.

The strength of the first derivative regularizers along spatial dimension is controlled by η s ≥ 0. The strength of the zero-attracting regularizer is controlled by η z ≥ 0.

The spectral regularization term D λ x k-q+1 2 promotes spectral smoothness which corresponds to the very nature of NIR hyperspectral images of biological material such as wood [START_REF] Tsuchikawa | A review of recent application of near infrared spectroscopy to wood science and technology[END_REF]. Matrix D λ is a first-order filtering operator along the spectral dimension weighted by the coefficients {c p } P -1 p=1 . It is defined as:

(3.6) D λ (diag(c 1 , . . . , c P -1 ) T P ) ⊗ I N .
The parameter η λ controls the strength of the spectral smoothness penalty term.

The spectral weights c p were introduced to provide some additional flexibility to the algorithm. Basically, the strength of the weight c p promotes the smoothness around the p-th spectral band. In our case, the spectral data are sampled on a regular grid.

hence, c 1 , . . . , c P -1 are all set to 1. however, in the case where the spectral data would be recorded on two non adjacent spectral bands, smoothness between these bands should not be promoted and the corresponding weight should be set to a small value (zero).

The criterion (3.1) developed in [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] introduces a regularization using second order derivatives along spatial and time dimensions. An inspection of criterion (3.4) may lead to the assertion that no regularization term is acting along the time dimension (time-regularization). However, we prove in Appendix B that, for the denoising case, LMS-based algorithm implicitly introduces a time regularization and presents a time delay in the estimation. A formal link between µ and the regularization strength is given. A small value µ results in strong regularization strength and large delay along the time dimension. Experimental results (Subsection 5.1) confirm that the interpretation is valid for the deconvolution case.

A subgradient of (3.4) is given by:

(3.7) ∇J (x k , . . . , x k-Q-L+2 ) ∂J ∂x k , . . . , ∂J ∂x k-Q+1 , 0 N ×1 , . . . , 0 N ×1
where 0 I×J denotes the I × J zero matrix. Zero terms in the subgradient indicate that past estimates xk-Q , . . . , xk-Q-L+2 are no longer updated. We now derive the sliding-block regularized LMS algorithm. Consider the vectorized data:

(3.8) x ′ k col{x k-q+1 } Q+L-1 q=1 , y ′ k col{y k-q+1 } Q+L-1 q=1
.

A valid subgradient for |x| is sign(x) [START_REF] Eksioglu | Sparsity regularised recursive least squares adaptive filtering[END_REF] where the sign function is defined as sign(x) = 0 for x = 0, and sign(x) = x/|x| otherwise. Approximating the subgradient in (3.7) by its instantaneous value yields:

(3.9) ∇J (x ′ k ) = -2Φ (y ′ k -Gx ′ k ) + η z Γsign(x ′ k ) + η s Λ ⊤ s sign(Λ s x ′ k ) + 2η λ Λ ⊤ λ Λ λ x ′ k Matrices Φ, G and Γ are of size (Q + L -1)P N × (Q + L -1)P N . Matrix Φ is defined by: (3.10) Φ Φ ⊤ 11 0 QP N ×(L-1)P N 0 (L-1)P N ×(Q+L-1)P N
where Φ 11 is the matrix of size QP N × QP N defined as:

Φ 11     H 1 • • • H Q . . . . . . 0 H 1     (3.11)
in which H ℓ = 0 P N ×P N for ℓ > L. It is necessary to define the matrix Φ 11 as above to properly account for the cases Q < L and Q ≥ L. Matrix G has the form:

(3.12) G       H 1 • • • H L 0 . . . . . . 0 H 1 • • • H L 0 (L-1)P N ×(Q+L-1)P N       Φ 11 G 12
0 (L-1)P N ×(Q+L-1)P N , (see Appendix A for more details). Matrix Γ is defined by:

(3.13) Γ I QP N 0 QP N ×(L-1)P N 0 (L-1)P N ×(Q+L-1)P N .
The first-order derivative filters for spatial and spectral dimensions are:

(3.14) Λ s [I Q ⊗ D s 0 QP (N -1)×(L-1)P N ] [Λ s,1 0 QP (N -1)×(L-1)P N ], (3.15) Λ λ [I Q ⊗ D λ 0 Q(P -1)N ×(L-1)P N ] [Λ λ,1 0 Q(P -1)N ×(L-1)P N ].
Finally, the SBR-LMS algorithm for hyperspectral image deconvolution is given by:

(3.16) x′ k+1 = Ωx ′ k +µΦ (y ′ k -Gx ′ k )-ρ z Γsign(x ′ k )-ρ s Λ ⊤ s sign(Λ s x′ k )-µη λ Λ ⊤ λ Λ λ x′ k
where ρ z = µη z /2, ρ s = µη s /2 and µ is a step size parameter that controls the trade off between convergence rate and algorithm stability. Matrix Ω is given by: (3.17)

Ω    I QP N 0 QP N ×(L-1)P N 0 I P N 0 0 P N ×P N 0 0 I (L-2)P N 0 (L-2)P N ×P N    I QP N 0 QP N ×(L-1)P N Ω 21 Ω 22 .
The upper part of matrix Ω corresponds to the set of updated variables, and the lower part allows to shift the past estimates. The final result xk-Q+2 is obtained by selecting the Q-th block of vector x′ k+1 , that is,

(3.18) xk-Q+2 = Sx ′ k+1
where S [0 P N ×(Q-1)P N , I P N , 0 P N ×(L-1)P N ]. It is worth to mention that the proposed algorithm is different from the standard block-LMS algorithm for which the output x k is updated only once for every block of size Q. On the contrary, in the proposed algorithm, to account for the causality issues discussed in Subsection 2.3,

x k is updated Q times.
Finally, depending on the hyperparameter values, different LMS-like algorithms can be defined as presented in Table 1.

Table 1: The different LMS-like algorithms and corresponding hyperparameters We also specify the particular case of using separable convolution kernel in form of

µ Q ηz ηs η λ ZA-LMS > 0 1 > 0 0 0 SB-LMS > 0 > 1 0 0 0 SBZA-LMS > 0 > 1 > 0 0 0 SBR-LMS > 0 > 1 > 0 > 0 > 0 3
H ⋆p = g p f p⊤ where f p = [f p L , . . . , f p 1 ] ⊤ and g p = [g p M , . . . , g p 1 ]
⊤ which yields lower computational burden. The validity of the separable kernel assumption is mainly controlled by the calibration of the optical systems. If the optical elements are misaligned, then the kernel will be non separable. Otherwise, a separable Gaussian kernel is a good and widely used model (see [START_REF] Goodman | Introduction to Fourier optics[END_REF] for details).

The computational complexity per iteration of SBR-LMS is given in Table 2,

where the costs of both direct and frequency domain implementations of the convolution are given. The approximate costs in Table 3 give the number of multiplications obtained by setting Q = L = M and retaining the dominant complexity. In the separable kernel case, both direct and frequency domain implementations yield lower complexity than for a non-separable kernel. In the direct domain, the complexity is linear w.r.t. the number of unknowns P N . For Q ≫ log 2 N , we can reduce the complexity for the non-separable kernel case by applying the algorithm in the frequency domain while for separable kernel case, both direct domain and frequency domain yield a computational cost of the same order. As expected, when Q ≪ log 2 N , the frequency domain implementation provides no benefit.

10

This manuscript is for review purposes only.

Ωx ′ k + µΦ y ′ k -Gx ′ k direct domain,
non-separable kernel

P Q L + Q + 1 2 M N - M -1 2 P QL (M -1)(N - M -1 2 ) + N + P Q Q + 1 2 (M -1)(N - M -1 2 ) + Q -1 2 N Ωx ′ k + µΦ y ′ k -Gx ′ k direct domain,
separable kernel

P Q L + Q + 1 2 N + 2P QM N - M -1 2 P Q L + Q -1 2 N + 2P Q(M -1)(N - M -1 2 ) Ωx ′ k + µΦ y ′ k -Gx ′ k frequency domain, non-separable kernel 4P Q L + Q + 1 2 N + 2P Q(L + 1)N log 2 N P Q(2L + Q -1)N + 2P Q(L + 1)N log 2 N Ωx ′ k + µΦ y ′ k -Gx ′ k frequency domain, separable kernel P Q(L + Q + 1)N + 2P QN (2log 2 N + 4) P Q(L + Q)N + 4P QN log 2 N ρzΓsign(x ′ k ) 0 0 ρsΛ ⊤ s sign(Λs x′ k ) 0 2(N -1)P Q µη λ Λ ⊤ λ Λ λ x′ k 4N (P -1)(Q + L -1) 2N P (P -1)(Q + L -1) 11
This manuscript is for review purposes only.

Algorithm 3.1 SBR-LMS algorithm for hyperspectral image deconvolution

Data: Y ∈ R N ×K×P , H ⋆p = [h p L , . . . , h p 1 ], ∀p = 1, . . . , P Result: X ∈ R N ×K×P Initialization Q, µ, ρ z , ρ s , η λ , c 1 , . . . , c P -1 and T N h p ℓ = 0 N ×1 , ∀ℓ = L + 1, . . . , Q, ∀p = 1, . . . , P ;
Selection of parameters X = 0 N ×K×P , c 0 = c P = 0 for k = L + Q -1 : K -1 do for q = 1 : Q do for p = 1 : P do Non-separable convolution kernel:

err p k-q+1 = y p k-q+1 - L ℓ=1 h p ℓ * x p k-q-ℓ+2 grad = q i=1 h p q-i+1 * err p k-i+1
Separable convolution kernel:

err p k-q+1 = y p k-q+1 -g p * L ℓ=1 f p ℓ x p k-q-ℓ+2 grad = g p * q i=1 f p q-i+1 err p k-i+1
where * is the 1D convolution operator. 2) can be written as:

x p k-q+2 =x p k-q+1 + µ grad -ρ z sign(x p k-q+1 ) -ρ s T T N sign(T N x p k-q+1 ) -µη λ (c 2 p-1 + c 2 p )x p k-q+1 -c 2 p-1 x p-1 k-q+1 -c 2 p x p+1 k-q+1 ; end end end
O(Q 3 P N ) O(Q 2 P N log 2 N ) Separable kernel O(Q 2 P N ) O(QP N (Q + log 2 N )) 4
y ′ k = Gx ′ k + z ′ k with z ′ k a
zero-mean measurement noise of covariance σ 2 z I (Q+L-1)N P . Taking the expectation of both sides of (3.16), the mean sliding block vector evolves according to the following recursion:

(4.1) E{x ′ k+1 } = AE{x ′ k } + µΦGx ′ k -ρ z Γ E{sign(x ′ k )} -ρ s Λ ⊤ s E{sign(Λ s x′ k )}
where A Ω -µΦG -µη λ Λ ⊤ λ Λ λ . The mean vector can be then obtained by using:

(4.2) E{x k-Q+2 } = SE{x ′ k+1 },
The main difficulty in (4.1) lies in evaluating the expectation of the sign function.

This point will be discussed later.

Consider the error:

(4.3) ǫ k = S(y ′ k -Gx ′ k ) = S (z ′ k -G(x ′ k -x ′ k )) ,
as z ′ k is assumed independent of other variables, the mean-squared error (MSE) can be expressed as:

E ǫ k 2 =N P σ 2 z + trace(G ⊤ S ⊤ SGE{x ′ k x′⊤ k }) + x ′ k ⊤ G ⊤ S ⊤ SGx ′ k -2x ′ k ⊤ G ⊤ S ⊤ SGE{x ′ k }.
It follows that E{x ′ k x′⊤ k } can be updated as:

E{x ′ k+1 x′⊤ k+1 } =AE{x ′ k x′⊤ k }A ⊤ + µ 2 σ 2 z ΦΦ ⊤ + µΦGx ′ k E{x ′⊤ k }A ⊤ + µAE{x ′ k }x ′⊤ k G ⊤ Φ ⊤ + µ 2 ΦGx ′ k x ′⊤ k G ⊤ Φ ⊤ -ρ z (U 1 + U ⊤ 1 + U 2 + U ⊤ 2 ) + ρ 2 z U 3 -ρ s (U 4 + U ⊤ 4 + U 5 + U ⊤ 5 ) + ρ 2 s U 6 + ρ s ρ z (U 7 + U ⊤ 7 )
where

U 1 = AE{x ′ k sign(x ′ k ) ⊤ }Γ ⊤ , U 2 = µΦGx ′ k E{sign(x ′ k ) ⊤ }Γ ⊤ , U 3 = ΓE{sign(x ′ k )sign(x ′ k ) ⊤ } Γ ⊤ , U 4 = AE{x ′ k sign(Λ s x′ k ) ⊤ }Λ s , U 5 = µΦGx ′ k E{sign(Λ s x′ k ) ⊤ }Λ s , U 6 = Λ ⊤ s E{sign(Λ s x′ k )sign(Λ s x′ k ) ⊤ }Λ s , U 7 = ΓE{sign(x ′ k )sign(Λ s x′ k ) ⊤ }Λ s .
Again, the main difficulty lies in the evaluation of the expectations in form of E{sign(u) sign(v)} and E{u sign(v)}. An analysis of the stochastic behavior of ZA-LMS is carried out in [START_REF] Chen | Transient performance analysis of zeroattracting lms[END_REF] where exact expressions of these expectations are derived under a Gaussian assumption on the error vectors. however, due to the size of the covariance matrix (4.4), we decided to resort to use the following approximations [START_REF] Zhang | Transient analysis of zero attracting NLMS algorithm without gaus-sian inputs assumption[END_REF]:

E{sign(u)} ≈ sign(E{u}), (4.4) E{u sign(v) ⊤ } ≈ E{u}sign(E{v}) ⊤ , (4.5) E{sign(u) sign(v) ⊤ } ≈ sign(E{u})sign (E{v}) ⊤ (4.6) 13
This manuscript is for review purposes only.

which lead to sufficiently accurate models. When |u| and |v| are large, the approximations are quite good meaning that the mean behavior is well approximated when the dynamic of the image to restore is high enough. When they are close to zero, the approximations are of lower quality. However, following [START_REF] Chen | Regularized least-mean-square algorithms[END_REF], these approximations yield a worst case transient behavior model. In addition, this approximation does not assume any statistical model on u and v. Only their statistical independence is required.

The mean and mean-squares transient models will be used to study the SBR-LMS behavior in both steady-state and non-stationary cases (Subsection 5.1). In the next section, to study the stability condition of the algorithm, we focus on the steady-state case.

Stability condition.

We shall now discuss the stability of the proposed algorithm in the mean and mean-squares sense. Consider a steady-state solution for which: (4.7)

x p k = xp , ∀ k.

The error-vector v p k is defined as: k can be updated as:

(4.8) v p k = xp -xp k with its vectorized version v k col{v p k } P p=1 (4.9) v ′ k col{v k-q+1 } Q+L-1 q=1 . (4.
(4.11) v ′ k+1 = Av ′ k -µΦz ′ k + (I -Ω + µη λ Λ λ ⊤ Λ λ )x ′ + ρ z Γsign(x ′ k ) + ρ s Λ s ⊤ sign(Λ s x′ k ) = Av ′ k -µΦz ′ k + µη λ Λ λ ⊤ Λ λ x′ + ρ z Γsign(x ′ k ) + ρ s Λ s ⊤ sign(Λ s x′ k )
where x col{x p } P p=1 and x′ [x ⊤ , . . . , x⊤ ] ⊤ of size (Q + L -1)P N × 1. We consider the following partition of the matrix A:

A = I 0 Ω 21 Ω 22 -µ Φ ⊤ 11 0 0 0 Φ 11 G 12 0 0 -µη λ Λ ⊤ λ,1 Λ λ,1 0 0 0 A 11 A 12 Ω 21 Ω 22 (4.12)
where Φ 11 is defined by (3.11), Ω 21 and Ω 22 are defined by (3.17),

A 11 = I - µ(Φ ⊤ 11 Φ 11 + η λ Λ ⊤ λ,1 Λ λ,1 ), A 12 = -µΦ ⊤ 11 G 12 .
Mean stability. We now consider the mean stability of (4.11). Taking the expectation of both sides, the mean-error vector can be expressed as:

(4.13) E{v ′ k+1 } = AE{v ′ k } + d k where (4.14) d k µη λ Λ ⊤ λ Λ λ x′ + ρ z E{sign(x ′ k )} + ρ s Λ ⊤ s E{sign(Λ s x′ k )}.
Let E{v ′ 1 } be the initial condition. Then we obtain:

(4.15) E{v ′ k+1 } = (A) k E{v ′ 1 } + k-1 n=0 (A) n d k-n
The Euclidean norm of d k is bounded. Indeed, we have:

(4.16) ||d k || ≤ µη λ ||Λ ⊤ λ Λ λ x′ || + ρ z ||E{sign(x ′ k )}|| + ρ s ||Λ ⊤ s E{sign(Λ s x′ k )}||.
The first term ||Λ ⊤ λ Λ λ x′ || corresponds to the initial condition and is bounded if x′ is bounded. The last two terms are also bounded:

||E{sign(x ′ k )}|| ≤ QP N (4.17) ||Λ ⊤ s E{sign(Λ s x′ k )}|| ≤ QP (4N -6). (4.18)
Following the arguments reported in [START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF], if the spectral norm of A is less than 1, the series ∞ n=0 (A) n d k-n is convergent. Thus, the convergence of (4.15) is entirely controlled by the spectral norm of A. Because matrix A is not symmetric, standard stability results cannot be applied straightforwardly. Let D be the closed unit disc in the complex plane. We say that the matrix A is discrete stable (see [START_REF] Buchevats | Linear Discrete-Time Systems[END_REF]) if δ(z) = det(I -zA) = 0, ∀z ∈ D. Using the Schur determinant formula, it can be shown that:

(4.19) δ(z) = det(I -zA 11 ) • det(I -zΩ 22 -z 2 Ω 21 (I -zA 11 ) -1 A 12 ).
Then, δ(z) = 0 if and only if the following two conditions hold:

det(I -zA 11 ) = 0, ∀z ∈ D (4.20) det(I -zΩ 22 -z 2 Ω 21 (I -zA 11 ) -1 A 12 ) = 0, ∀z ∈ D (4.21)
The condition in (4.20) is fulfilled if all the eigenvalues of the symmetric matrix A 11 lie inside the unit disc, i.e., its spectral radius r(A 11 ) < 1. The matrix Φ ⊤ 11 Φ 11 + η λ Λ ⊤ λ,1 Λ λ,1 being the sum of a positive semi-definite matrix and a positive definite matrix, it is positive definite; thus, standard stability result applies from which we conclude that A 11 is stable iff:

(4.22) 0 < µ < 2 r(Φ ⊤ 11 Φ 11 + η λ Λ ⊤ λ,1 Λ λ,1 )
.

In practice, it is not easy to check the second condition in (4.21) for all values of z in the unit disc. However, in the case of Gaussian convolution kernels, we checked experimentally that the stability of A depends only on that of A 11 . For example, Figure 4 displays the spectral radius of A as a function of the tuning parameters µ and η λ for two Gaussian kernels. The red curve is obtained using (4.22) which corresponds to the stability limit of A 11 . The region where the spectral radius of A is less than 1 (i.e., A is stable) is highlighted in green. One can observe that for different Gaussian kernels, except for some numerical errors, the stability region of A is the same as that of A 11 . In conclusion, we conjecture that for Gaussian filters, the stability condition of A is also given by (4.22). Numerical experiments seem to indicate that the stability of A is also controlled by the stability of A 11 for a larger class of low-pass filters, including hanning, hamming and rectangular windows. We also found counter examples for which the stability of A 11 does not guarantee the stability of A. It is worth mentioning that the stability region depends on both µ and η λ and that increasing η λ necessitates to lower µ to guarantee the stability of the algorithm. Mean-squares stability. We now address the stability in the mean-squares sense.

First, define the zero-mean misalignment vector w k as:

(4.23) w k = v ′ k -E{v ′ k },
then using (4.11) and (4.13), w k can be updated as: The covariance matrix of w k is defined as:

(4.27) W k = E{w k w ⊤ k }
Since E{w k } = 0, E{e zk } = 0 and E{e sk } = 0, W k can be updated as follows:

W k+1 =AW k A ⊤ + µ 2 σ 2 z ΦΦ ⊤ + ρ z AE{w k e z ⊤ k } + ρ z E{e zk w ⊤ k }A ⊤ + ρ 2 z E{e zk e z ⊤ k } + ρ s AE{w k e s ⊤ k } + ρ s E{e sk w ⊤ k }A ⊤ + ρ 2 s E{e sk e s ⊤ k } + ρ s ρ z E{e sk e z ⊤ k } + ρ s ρ z E{e zk e s ⊤ k } (4.28)
From (4.8) and (4.23) we have:

E{w k e z ⊤ k } = E{x ′ k }E{e z ⊤ k } -E{x ′ k e z ⊤ k } = -E{x ′ k e z ⊤ k }, (4.29) E{w k e s ⊤ k } = E{x ′ k }E{e s ⊤ k } -E{x ′ k e s ⊤ k } = -E{x ′ k e s ⊤ k }. (4.30)
Thus, we observe:

trace(W k+1 ) =trace(A ⊤ AW k ) + µ 2 σ 2 z trace(ΦΦ ⊤ ) + 2ρ s ρ z E{e s ⊤ k e zk } -2ρ z trace(AE{x ′ ⊤ k e zk }) + ρ 2 z E{e z ⊤ k e zk } -2ρ s trace(AE{x ′ ⊤ k e sk }) + ρ 2 s E{e s ⊤ k e sk }. (4.31) 16
This manuscript is for review purposes only.

It is straightforward to verify that e zk and e sk are bounded. With

E{x ′ k } = x′ - E{v ′ k }, the terms E{x ′ ⊤ k e zk } and E{x ′ ⊤ k e sk } converge if E{v ′ k } converges which is
true if A is stable and x′ is bounded. Thus to have the stability of (4.28), only the stability of its first term on the right hand side has to be studied (see for example [START_REF] Shi | Convergence analysis of sparse lms algorithms with l1-norm penalty based on white input signal[END_REF]). By using the same argument as in [START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF], the mean-squares stability of (4.28) is guaranteed if the mean stability holds. The next experiment intends to show that the proposed model is still valid in the non-stationary case. The results are reported in Figure 6. The n-th pixel at time k of the 2D simulated image is defined as:

(5.2) and Figure 8(d) (with spectral regularization), it appears that spectral regularization can help to recover objects whose amplitude is very small (at a given wavelength) as pattern 1 in the second band. However, the spectral regularization may also induce artifacts: for example, pattern 2 appears in the second band when the spectral regularization is used while it is not present in the original image. This manuscript is for review purposes only. The convolution filter H ⋆ was assumed to be a low-pass Gaussian filter of size (9 × 9).

x (nonstat.) nk = x (stat.) n + |x (stat.) n | 10 sin 2π T o k + 2π n -1 K with T o =
It is invariant with respect to spectral band p. A Gaussian noise was then added to the blurred image. The noise level was the same for all bands.

The SBR-LMS results are compared to an off-line deconvolution algorithm proposed in [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] which corresponds to the Tikhonov approach with ℓ 2 spatial and spectral regularizers. It is also compared to the BT approach presented in Subsection 3.1 whose block size is the same as the block size of the SBR-LMS. Note that better performance would certainly be achieved by using much more refined methods enforcing additional constraints (non-negativity, edge preserving) at the price of an increased computational burden. That is why we decided to restrict the comparison to deconvolution approaches having a computational burden similar to that of the SBR-LMS.

Three versions of SBR-LMS algorithms are considered for performance evaluation: SBR-LMS without spatial and zero-attracting regularizations, SBR-LMS with spatial 23 This manuscript is for review purposes only.

regularization but without zero-attracting regularization, SBR-LMS with spatial and zero-attracting regularizations.

As a performance measure, we consider the MSE as a function of the SNR and the goal is to assess the performance of the proposed SBR-LMS and to compare them to those of alternative approaches. As shown in the previous section, the optimal hyperparameter values are depending on the noise level. Thus, comparing the performances of different methods would require to determine for each noise level the set of optimal parameters by an exhaustive grid search which, due to the large number of hyperparameters (2 for Tikhonov, 3 for Block Tikonov and up to 5 for SBR-LMS), is out of reach in reasonable time. To handle this problem, while keeping a fair comparison between the different methods, we propose the following performance evaluation: the strength of the spectral regularization is fixed to a small value (the same for all methods) and the block size Q is fixed to 11 (BT, SBR-LMS). All the other parameters are fixed by a grid search yielding the minimum MSE for a fixed SNR= 12 dB. These parameter values are then used for all other noise levels. Figure 12 shows that SBR-LMS without spatial and zero-attracting regularization terms yields the highest MSE for all noise level. This can be attributed to the fact that, in that case, only a time regularization is used. For high SNR, Tikhonov and BT give the lowest MSEs. however, when SNR is lower than 12 dB, SBR-LMS with spatial regularization has almost the same performance as the Tikhonov approach and better performance than BT algorithm. Adding a zero-attracting regularization improves the performance at low SNR even more. In conclusion, if we focus on sequential deconvolution, BT has to be chosen if the SNR is high while SBR-LMS should be preferred for low SNR. The spectral response of the conveyor (background) was estimated from data in an area of size 140×140×28. It was then subtracted from each pixel of the hyperspectral image. The imaged objects put on the conveyor are pieces of woods (wood wastes).

The convolution filter was estimated from data (using a calibration target) to be a Gaussian kernel of size 7 × 7 with FWHM of 3 points. No dependence of the convolution kernel w.r.t. the wavelength was observed.

First, we compare different algorithms in the case of an image acquired using a large integration time (2.146 ms) resulting in a high SNR. This corresponds to the situation where Tikhonov and BT are expected to yield the best results. Figure 13 presents the deconvolution result obtained on the whole real hyperspectral image (only 3 slices corresponding to wavelengths 1058.64 nm, 1198.12 nm and 1479.78 nm are shown). Figure 14 The effective implementation of such method is hampered by the choice of the regularization parameters. In general, this choice is made by successive trials which can be highly time consuming. Future works will focus on the automatic learning of hyperparameters. A first idea could be to learn the optimal parameters in a way similar to what is done in [START_REF] Song | Regularization parameter estimation for non-negative hyperspectral image deconvolution[END_REF] on a typical sample and to use these fixed parameters for the online processing. A more promising approach would be to perform the online learning of these parameters: for example, [START_REF] Jin | Adaptive parameters adjustment for group reweighted zero-attracting[END_REF], a variable-parameter algorithm is proposed to adjust optimally the algorithm parameters of ZA-LMS. This will avoid the necessity of having representative image samples and will meet the requirements

25
This manuscript is for review purposes only. For Q ≥ L, matrix G can be partitioned as follows:

730 (A.1) G             H 1 • • • H L 0 . . . . . . H 1 • • • H L . . . . . . H 1 • • • H L 0 (L-1)P N ×(Q+L-1)P N             , 731 27
This manuscript is for review purposes only. and for Q < L, we have the form:

732 (A.2) G       H 1 • • • H Q • • • H L 0 . . . . . . H 1 • • • H Q • • • H L 0 (L-1)P N ×(Q+L-1)P N      
.

733 Thus, we can conclude that for both cases Q ≥ L and Q < L, matrix G can be 734

Fig. 1 :

 1 Fig. 1: Data acquisition by hyperspectral imaging systems

Fig. 2 :

 2 Fig. 2: Structure of a hyperspectral data cube Y

Fig. 3 :

 3 Fig. 3: Data acquisition with (a) low and (b) high spatial sampling rate. high spatial sampling rate results in a blurring of the image.

3. 2 .

 2 Sliding-block regularized LMS (SBR-LMS). Consider the problem of sequentially estimating x k-Q+1 based on the delayed observations (ỹ k , . . . , ỹk-Q+1 ) and past estimates xk-Q , . . . , xk-Q-L+2 . Due to the non-causality of the estimation process, we propose to determine the future estimates with respect to time k

. 3 .

 3 Algorithm implementation and computational cost. Algorithm 3.1 details the SBR-LMS implementation. It should be noted that the storage of matrix H p ℓ is avoided and that only 1D convolutions are required. Two different implementations of the convolution are possible: direct (spatial) domain or frequency domain.

  . Transient behavior analysis. In this section, we derive the transient behavior model of the deconvolution algorithm (3.16) both in the mean and mean-squares sense. The analysis provides us with an important guidance for selecting the stepsize parameter, and also allows us to study the influence of the other parameters (Subsection 5.2). The stability condition of the algorithm is also studied.

4. 1 .

 1 Mean and mean-squares transient behavior model. Using the notations of Subsection 3.2, Equation (2.

10 )

 10 By combining (3.16), (4.8) and (4.10), v ′

  L = 11, Q = 12

Fig. 4 :

 4 Fig. 4: Spectral radius of A (in log scale) as a function of µ and η λ . The green area corresponds to the values of µ and η λ for which A is stable. The red curve corresponds to the stability limit of A 11 .

(4. 24 )

 24 w k+1 = Aw k -µΦz ′ k + ρ z e zk + ρ s e sk where e zk = sign(x ′ k ) -E{sign(x ′ k )}, (4.25) e sk = Λ ⊤ s (sign(Λ s x′ k ) -E{sign(Λ s x′ k )}). (4.26)

5. Experimental results. 5 . 1 . 2 F 2 F.

 5122 Figure 5(b).

  100. The convolution filter and noise level are the same as before. Q is set to 5, ρ z and ρ s are set to 0.2 • 10 -3 . The simulated images (in black dots) have to be compared with the superimposed experimental (blue) and theoretical (red) curves. Figure6(a), Figure 6(b) and Figure 6(c), Figure 6(d) confirm that the step size influences the convergence speed and the estimation variance. A large value for µ results in an increase of the convergence speed and estimation variance. Also, the delay mentioned in Subsection 3.2 and Appendix B can be observed on these figures. Increasing the step size (Figure 6(c)) results in a decreasing of the delay. 5.2. Effects of the parameters. To assess the influence of the parameters Q, µ, ρ s and ρ z , numerical simulations were conducted on a 2D image. The original image shown in Figure 7(a) was blurred by a Gaussian kernel of size 15 × 15 with Evolution of E{x k },ρs = ρt = ρz = 0 Evolution of MSE, ρs = ρz = 0 Evolution of E{x k },ρz = 1.5 • 10 -3 ,ρs = 0 Evolution of MSE,ρz = 1.5 • 10 -3 , ρs = 0 Evolution of E{x k },ρz = 0, ρs = 1.5 • 10 -3 Evolution of MSE,ρz = 0, ρs = 1.5 • 10 -3

Fig. 5 :

 5 Fig. 5: Transient behavior model in stationary environment

Figure

1 Fig. 6 :

 16 Fig. 6: Transient behavior model validation in non-stationary environment

Figure 8

 8 Figure 8 shows the influence of the spectral regularization parameter η λ on a simulated hyperspectral image. Original and blurred noisy images corresponding to three different wavelengths are shown in Figure 8(a) and Figure 8(b), respectively. The blurring kernel and the noise level are the same as those in Figure 7(b). By comparing he estimated results shown in Figure 8(c) (without spectral regularization)

  We can assess the effect of the parameter Q on the convergence of the SBR-LMS by examining the transient behavior of the MSE. The simulated image is the same as the one in Subsection 5.1 (2D image whose lines are constant over time k). It was blurred by Gaussian filter of size 5 × 5 with FWHM equal to 3 and the SNR of the 19 This manuscript is for review purposes only.

Fig. 7 :

 7 Fig. 7: Influence of different hyperparameters

Fig. 8 :

 8 Fig. 8: Influence of spectral regularization parameter

Fig. 9 :Fig. 10 :

 910 Fig. 9: Evolution fo MSE with different values of Q

  Fig. 11: MSE as function of k for different values of ρ z

3 Tikhonov, ηs = 8 ,Fig. 12 :

 3812 Fig. 12: Performances of the hyperspectral image deconvolution algorithms

  Figure 14(b). Figure13(c) and Figure14(c) correspond to the Tikhonov approach with ℓ 2 spatial and spectral regularizers. Figure13(d) and Figure14(d) correspondto the BT approach. The hyperparameters of the Tikhonov approaches were estimated by the minimum distance criterion (MDC) proposed in[START_REF] Song | Regularization parameter estimation for non-negative hyperspectral image deconvolution[END_REF]. The FWHB of the kernel being equal to 3 points, the blurring of the real images is quite limited. However, an increased resolution allowing to recover the fine texture of the objects can be observed on the zoomed view restored images in Figure14(b)-Figure14(d). More importantly, Figures Figure13and Figure14show that at high SNR, the SBR-LMS yields similar

  Figures Figure13and Figure14show that at high SNR, the SBR-LMS yields similar performances as those of the reference methods (Tikhonov and BT).The second image in Figure15(a) corresponds to the same observed scene with an integration time of 0.013 ms which results in a low SNR. The image restored with SBR-LMS (µ = 0.5, Q = 9, ρ z = 0.1, ρ s = 0.05, η λ = 0.001) is shown in Figure 15(b). Results obtained by Tikhonov and BT methods, whose hyperparameters are estimated by the MDC approach, are shown in Figure 15(b) and Figure 15(c), respectively. The bias introduced by the regularization terms can be observed in Figure 15(b) since the dynamic range is lower than that of Tikhonov methods. It clearly appears that the noise level on the background estimated by SBR-LMS is much lower than that of the Tikhonov approaches. Due to the noise level, all approaches cannot reveal fine structures accurately, but the edges are well preserved. 6. Conclusions. In this work, the online deconvolution problem of hyperspectral images collected by industrial pushbroom imaging systems was addressed. The contribution of this work is the derivation of the SBR-LMS algorithm that allows the fast slice-by-slice hyperspectral image restoration, accounting for convolution kernel non-causality and low SNR issues. The key feature of the proposed approach relies on the possibility to include non-quadratic regularization terms while maintaining a computational complexity linear w.r.t. the number of unknowns. The transient behavior model of the algorithm was analyzed; it allows to assess the influence of each regularization parameter. Experimental results on both simulated and real hyperspectral images showed that the SBR-LMS algorithm has good noise removal and deblurring capabilities, especially at low SNR which is the relevant case for industrial imaging systems.

Fig. 13 :Fig. 14 :

 1314 Fig. 13: hyperspectral image restoration at 3 wavelengths.

Fig. 15 :

 15 Fig. 15: hyperspectral image restoration at 3 wavelengths.

  ×K is the blurred spatial image for each wavelength λ p , X p ∈ R N ×K is the image to restore, H ⋆p ∈ R M×L is a convolution kernel (filter), and Z p ∈ R N ×K

	is an additive i.i.d. noise. The symbol * stands for the 2D convolution operator. First we derive a sequential causal formulation of model (2.1). Without loss of generality,
	we focus on the sequential model for 2D images. The image Y p , collected online, can be represented as a sequence of vectors y p k := [y p 1,k , . . . , y p N,k ]

⊤ , with k = 1, . . . , K,

Table 2 :

 2 Computational cost per iteration of recursion(3.16)

	Operation	Multiplications	Additions

Table 3 :

 3 Approximate computational cost per iteration for the convolution in(3.16)

	Direct domain	Frequency domain
	Non-separable kernel

For simplicity, L is assumed to be odd.

For simplicity, M is assumed to be odd.

For simplicity, Q is chosen to be odd.

This manuscript is for review purposes only.

This work has been supported by the FUI AAP 2015 Trispirabois Project, the Conseil Régional de Lorraine, the GDR ISIS CNRS, and the CNRS Imag'in ALOhA Project. The work of J. Chen was supported in part by National Natural Science Foundation of China (NSFC) under grant 61671382 and in part by NSF of Shenzhen under grant JCYJ2017030155315873.

† Centre de Recherche en Automatique de Nancy (CRAN),

in [START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF] could be considered. For different types of materials and/or spectroscopies, 722 other spectral regularizations could be used such as TV-like spectral regularization 723 proposed in [START_REF] Aggarwal | Hyperspectral image denoising using spatio-spectral total variation[END_REF] for promoting piecewise constant spectra reconstruction. Finally, 724 accounting for the low-rank structure of the data to restore (which results from a 725 non-negative linear mixing model of the data to restore), a joint online deconvolution 726 and unmixing algorithm is worth being studied. This is expected to yield a very low 727 computational burden and accurate image restoration approach.

Appendix A. Structure of matrix G for both cases Q ≥ L and Q < L. Appendix B. SB-LMS interpreted as a time-regularized Tikhonov-like estimator. This appendix reveals that there is a hidden time-regularization introduced by the SB-LMS (that is the SBR-LMS where ρ s = ρ z = η λ = 0) whose strength is controlled by the step size µ. Setting the spatial and spectral regularization terms to zero in (3.16), the SB-LMS is defined as:

Denoting X′ (z) as the z-transform of x′ k , Equation (B.1) can be expressed as:

In another word, X′ (z) is the output of a linear filter:

with input Y ′ (z). It is well known that Tikhonov estimators result in a linear filtering of the data. Thus, the question we are addressing is to link these two linear filters.

To simplify the analysis, we consider the denoising problem for which the convolution kernel is set to be identity. In that case, the block length Q is set to 1 (which results in x′ k = xk ) and (B.2) can be simplified as:

Assuming a unit sampling step size, the forward transform is given by s = z -1 where s is the Laplace parameter. The Fourier transform on the filter F (z) = µ z -1 + µ is obtained by setting s = jω where ω is the angular velocity:

Combining the first-order approximation of arctan(ω/µ) ≈ ω/µ together with the series expansion of √ 1 + x 2 can be used to give a low-pass approximation of the filter: 2 . Following [START_REF] Moussaoui | Regularization aspects in continuous-time model identification[END_REF], any filter of the form (B.6) results from to the minimization of a criterion:

where D i is the ith-order (discrete) derivative matrix along the time dimension and k 0 = 1/µ. This means the LMS-based algorithms can be interpreted as delayed [START_REF] Hunt | Karhunen-Loeve multispectral image restoration, part I: Theory[END_REF] This manuscript is for review purposes only.

Tikhonov-like algorithms. The first point to mention is that LMS-based estimators present a time delay in the estimation which is (approximately) proportional to 1/µ. The smaller µ is, the larger delay is. This has to be opposed to Tikhonov estimators which are null-phase filters and do not introduce any delay. The second point is related to the regularization parameters β i (µ) which is proportional to 1/(2µ) 2i . A small value µ results in strong regularization strength along the time dimension.