G Tochon

L Drumetz

M A Veganzones

M Dalla Mura

J Chanussot

FROM LOCAL TO GLOBAL UNMIXING OF HYPERSPECTRAL IMAGES TO REVEAL SPECTRAL VARIABILITY

Keywords: Local spectral unmixing, extended linear mixing model, hyperspectral imagery

The linear mixing model is widely assumed when unmixing hyperspectral images, but it cannot account for endmembers spectral variability. Thus, several workarounds have arisen in the hyperspectral unmixing literature, such as the extended linear mixing model (ELMM), which authorizes endmembers to vary pixelwise according to scaling factors, or local spectral unmixing (LSU) where the unmixing process is conducted locally within the image. In the latter case however, results are difficult to interpret at the whole image scale. In this work, we propose to analyze the local results of LSU within the ELMM framework, and show that it not only allows to reconstruct global endmembers and fractional abundances from the local ones, but it also gives access to the scaling factors advocated by the ELMM. Results obtained on a real hyperspectral image confirm the soundness of the proposed methodology.

INTRODUCTION

Given a hyperspectral image (HSI), spectral unmixing is concerned with the identification of a set of signatures corresponding to pure macroscopic constituents, called endmembers and their respective proportions, termed fractional abundances, for all pixel spectra in the HSI. The linear mixing model (LMM) is a widely accepted assumption in which each pixel spectrum can be described as a linear combination of those endmembers and associated abundances [START_REF] Keshava | Spectral unmixing[END_REF]. The relative simplicity of the LMM has spawned numerous spectral unmixing techniques [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. However, in spite of its popularity, the LMM suffers from several drawbacks. In real life scenarios, the mixing process cannot always be modeled as linear (due to multiple scattering or intimate mixing effects). Thus, nonlinear unmixing methods have been proposed (see [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF] for a review). Moreover, the LMM implicitly assumes that each endmember can be properly represented by a single spectral signature, which is seldom true in practice due to the topography of the imaged scene and other physical aspects, including variable illumination conditions and the intrinsic variability of the materials (e.g., the amount of chlorophyll in vegetation). This endmember variability has received an increasing attention lately, and two main strategies have arisen in the literature [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF]. The first one models the endmember variability through some particular multivariate probability distribution [START_REF] Du | Spatial and spectral unmixing using the beta compositional model[END_REF]. The second approach typically describes the spectral variability of a given endmember as a set of different instances of this endmember signature. Spectral bundles [START_REF] Somers | Automated extraction of image-based endmember bundles for improved spectral unmixing[END_REF] are a typical example of such strategy where the spectral variability is not explicitly modeled. Lately, an extended linear mixing model (ELMM) has been proposed to address endmember variability while preserving the LMM framework [START_REF] Veganzones | A new extended linear mixing model to address spectral variability[END_REF], in which endmembers are allowed to vary from one pixel to another according to scaling factors. Local spectral unmixing (LSU) is a potential solution to alleviate this spectral variability phenomenon. As a matter of fact, physical causes inducing spectral variability can be assumed to be sufficiently moderate locally in the image for the classical LMM to hold. In [START_REF] Canham | Spatially adaptive hyperspectral unmixing[END_REF] for instance, the original image is first divided into several overlapping tiles, within which the spectral unmixing is locally conducted. To avoid blocking effects, another strategy was proposed in [START_REF] Drumetz | Binary partition tree-based local spectral unmixing[END_REF], where the LSU is performed on regions of a particular partition of the HSI. However, the main limitation of LSU techniques is that the induced endmembers and abundances are only defined locally in the image, and must somehow be post-processed to be interpretable at the whole image scale. In this work, we propose to combine the LSU approach with the ELMM framework and show that it allows to derive global endmembers and abundances from their local counterparts, while providing additional information related to the spectral variability. The rest of this paper is organized as follow: section 2 presents more in details the ELMM framework [START_REF] Veganzones | A new extended linear mixing model to address spectral variability[END_REF] and a LSU approach relying on a hierarchical representation of the HSI [START_REF] Drumetz | Binary partition tree-based local spectral unmixing[END_REF]. Section 3 describes the proposed methodology to connect the previous two strategies in order to have access to global endmembers, abundances, as well as spectral variability. Conducted experiments and results are presented in section 4, while conclusion and perspectives are drawn in section 5.

SPECTRAL UNMIXING AND VARIABILITY

Extended Linear Mixing Model

Let X ∈ R N λ ×Np be a HSI composed of Np pixels and comprising N λ spectral channels. Under the LMM assumption, each pixel spectrum x ∈ R N λ can be written as a linear combination of d endmember signatures ei, i ∈ {1, . . . , d} weighted by their corresponding fractional abundances φx,i, plus some additive noise η:

x = d i=1 φx,iei + η , (1)
with the fractional abundances subject to positivity (φx,i ≥ 0 ∀i) and sum-to-one ( d i=1 φx,i = 1) constraints. The geometrical interpretation of such model is that all data points x lie in a (d -1)-dimensional simplex whose d vertices are the endmembers ei. In the ELMM model [START_REF] Veganzones | A new extended linear mixing model to address spectral variability[END_REF][START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF] on the other hand, the endmember signatures ei are no longer constrained to be fixed, but can vary from one pixel to another according to some non negative scaling factor λi with respect to some reference signature e0 i . Following, equation (1) rewrites:

x = d i=1 φx,iλie0 i + η . (2)
As described in [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF], equation [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] states that each pixel lies on a simplex given by a cross-section of the convex hypercone spanned by the reference endmembers e0 i , and those parameters are estimated by solving an optimization procedure for all pixels of the HSI. The physical rationale behind the ELMM is to model variability due to illumination effects, but it also allows to recover information related to intrinsic variability [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF].

In any case, the quality of the unmixing process for a given pixel is frequently evaluated by the root mean square error (RMSE) between the true pixel spectrum x and the one reconstructed using the induced endmembers and abundances x:

(x, x) = 1 √ N λ x -x 2 .
(3)

Local Spectral Unmixing

In the LSU strategy, the spatial support sp(X) of the HSI X is first divided in a set of regions {R ⊆ sp(X)} which may overlap or not, but such that the LMM can be assumed to hold in each region R.

The unmixing of a pixel x belonging to a particular region R is then performed using only the information contained in this region:

x = d R i=1 φ R x,i e R i + η , (4)
where the dR endmembers e R i and their associated fractional abundances φ R

x,i have been induced using only the spectral information contained in region R. The major challenge of LSU approaches concerns the definition of a suited set of regions {R}. In this paper, we follow the work presented in [START_REF] Drumetz | Binary partition tree-based local spectral unmixing[END_REF] that makes use of a partition of the spatial support sp(X), π = {Ri ⊆ sp(X)} such that sp(X) = i Ri and Ri ∩ R j =i = ∅ extracted from a hierarchical representation of this HSI. More specifically, the HSI is first decomposed in a set of regions organized hierarchically (such that any two regions of this set are either disjoint or nested). The main advantage of hierarchical image decompositions is that they define regions at various scales ranging from very small (close to the pixels) to very large (up to the whole image spatial support) regions. The selection of regions with ideal scale to compose the resulting partition can then be formulated as an optimization problem conducted on the hierarchical structure [START_REF] Tarabalka | A marker-based approach for the automated selection of a single segmentation from a hierarchical set of image segmentations[END_REF]. Here, we exploit the idea presented in [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF], where LSU is first performed over all regions of the hierarchical representation, and its performance is assessed by the local reconstruction errors of the regions. The final obtained partition of the HSI is such that its overall reconstruction error is minimal with respect to all the other possible partitions contained in the hierarchical representation.

PROPOSED METHODOLOGY

As explained in the previous section 2.2, the starting point of the proposed methodology is a partition π = {Ri} of the spatial support of the HSI, where LSU procedures have been conducted on each region Ri of the partition. Following [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF], those regions have been designed to have minimal reconstruction errors. Thus,

dR i local endmembers e R i 1 , . . . , e R i d R i
are available in each region Ri, as well as their as-

sociated local fractional abundances φ R i x = φ R i x,1 , . . . , φ R i x,d R i
for all pixel spectra x belonging to Ri.

Endmember clustering

In a first step, all local endmembers are pooled together in a common set Eπ = {e R∈π i } dπ i=1 with dπ = R∈π dR being the total number of local endmembers that have been generated by the LSU approach. Then, this set Eπ is clustered into K clusters C1, . . . CK by means of some clustering algorithm, with cluster

C k = e R k,1 , . . . e R k,d C k
being composed of dC k local endmembers, originating from various regions R of the partition π. Finally, the centroid c k of each cluster C k is retrieved

c k = 1 dC k d C k i=1 e R k,i , (5)
and defined as the global endmember representing the k th cluster C k .

As opposed to all local endmembers e k,i belonging to this cluster, the centroid c k is expected to properly describe the macroscopic material associated with cluster C k across the whole image.

Global abundance retrieval

Once the global endmembers c1, . . . , cK have been defined, their associated global fractional abundances must be retrieved for all pixels spectra of the HSI. In particular, let x be such a pixel spectrum contained in region R of partition π, and let φ R x,1 , . . . , φ R

x,d R and [ψx,1, . . . , ψx,K ] be its local and global fractional abundances, respectively. Then, three possible cases may occur:

-No local endmember e R i has been clustered in C k . Thus, the k th macroscopic material represented by C k is not contained in x and ψ x,k = 0.

-There is a single local endmember e R i 0 belonging to cluster C k . Therefore, the proportion of this material shall not change within x, hence ψ x,k = φ R

x,i 0 . -There are several local endmembers e R in , n = 1, . . . , p grouped in the same cluster C k . In such situation, the material is locally variable within the region R (healthy and burnt grass for instance), but all contributions sum up with respect to the global instance of the material (being grass in the previous example). Thus, ψ x,k = p n=1 φ R x,in . All previous cases can be summarized as follows:

ψ x,k = d R i=1 1 {e R i ∈C k } φ R x,i , (6)
where

1 {e R i ∈C k } = 1 if e R i ∈
C k and 0 otherwise. Doing so for all pixels x of the HSI allows to reconstruct global fractional abundance maps.

Estimation of spectral variability

In the ELMM framework, each endmember ei is authorized to vary pixelwise with respect to some reference endmember e0 i according to some local scaling factor λi, as described by equation [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. Here, we take advantage of this idea by considering cluster centroids c k to be those reference endmembers and modeling all local endmembers belonging to this cluster C k as some scaled versions of c k :

e R i ∈ C k ⇒ e R i = λic k . (7)
The scaling factor λi associated with the local endmember e R i can be recovered in practice by least square regression between e R i and the centroid c k of cluster C k it belongs to:

λi = c T k c k -1 c T k e R i , (8)
with (•) T being the transpose operator. Besides, equation ( 8) guarantees the local scaling factor λi to be non negative.

Plugging equation ( 7) into equation ( 4) yields

x = d R i=1 φ R x,i λic k i + η , (9)
where ki ∈ {1, . . . , K} is the index of cluster e R i belongs to. Equation ( 9) can be rewritten as

x = K k=1 d R i=1 φ R x,i λi1 {e R i ∈C k } c k + η . ( 10
)
On the other hand, pixel x can also be decomposed with respect to the global ELMM framework as described by equation ( 2), namely

x = K k=1 ψ x,k ρ k c k + η = K k=1 d R i=1 φ R x,i 1 {e R i ∈C k } ρ k c k + η , (11)
with ρ k being the global scaling factor associated with centroid c k . Hence, it is possible to estimate the global scaling factor ρ k for pixel x as a weighted average of its local scaling factors λi and local abundances φ R x,i :

ρ k = d R i=1 φ R x,i λi1 {e R i ∈C k } d R i=1 φ R x,i 1 {e R i ∈C k } . ( 12
)
Note that, in the case where there is a single local endmember e R i 0 belonging to cluster C k , then ρ k = λi 0 is constant for all pixels x of region R. As a matter of fact, all pixels in R appear spatially homogeneous with respect to the material represented by C k . If there are at least two local endmembers belong to the same cluster C k on the other hand, then the global scaling factor ρ k varies pixelwise. Finally, in the case where

ψ x,k = d R i=1 φ R x,i 1 {e R i ∈C k } = 0 (that is, if pixel x does not contain the k th material), then ρ k is set to 1.

EXPERIMENTS

Experimental methodology

We apply the proposed methodology to the HSI acquired over the campus of the University of Houston in 2012 and described in [START_REF] Debes | Hyperspectral and LiDAR data fusion: Outcome of the 2013 GRSS data fusion contest[END_REF]. This image is composed of 340 × 320 pixels in spatial dimension, and comprises 144 bands. The study site features an urban area with a stadium, buildings, parking lots and roads, and some portions of grass and trees. A color composition of this HSI is presented by figure 1a. The partition π, input of the proposed methodology, is obtained following the procedure described in [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF]. First a spatial pre-processing of the HSI is conducted in order to mitigate the effects of potential outliers [START_REF] Zortea | Spatial preprocessing for endmember extraction[END_REF]. Then, a binary partition tree (BPT) representation [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF] is constructed using the spectral region model proposed in [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF] (modeling all regions by their local endmembers) and the endmember-based distance as merging criterion [START_REF] Grana | An endmember-based distance for content based hyperspectral image retrieval[END_REF]. A LSU procedure is conducted over each region of the BPT. The local intrinsic dimensionality dR of each region is estimated using random matrix theory [START_REF] Cawse-Nicholson | Determining the intrinsic dimension of a hyperspectral image using random matrix theory[END_REF], as it proved to be reliable when working over small regions [START_REF] Drumetz | Hyperspectral local intrinsic dimensionality[END_REF]. Local endmembers are estimated using the vertex component analysis (VCA) algorithm [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] and their fractional abundances are retrieved by non-negative sparse regression using a fractional semi-norm [START_REF] Meyer | Hyperspectral unmixing with material variability using social sparsity[END_REF]. The partition π extracted from the BPT structure is displayed by figure 1b and composed of 396 regions. It achieves a trade-off between low region-wise maximal reconstruction errors (penalizing large regions with potential high reconstruction errors, which may be caused by the invalidity of the LMM within the region) and simplicity (penalizing partitions with too many regions) (see [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF] for practical details).

All generated local endmembers are grouped in the set Eπ, eventually composed of 2957 individuals. This set Eπ is divided into K = 12 clusters following a multivariate Gaussian mixing model hypothesis, by application of the Expectation-Maximization algorithm [START_REF] Tarabalka | Spectralspatial classification of hyperspectral imagery based on partitional clustering techniques[END_REF]. Note that the total number of clusters K has been set empirically. According to the proposed methodology, all cluster centroids c k , k = 1, . . . , K are defined to be the ELMM reference endmembers. The associated global abundances ψ x,k and global scaling factors ρ k are retrieved for all pixel spectra x in the HSI following the procedures exposed in section 3.2 and section 3.3, respectively. In order to evaluate the performance of the proposed methodology, we also unmix the image following the classical LMM scenario: K "classical" endmembers and associated fractional abundances are globally induced over the image, using the same set-up as the LSU (namely VCA for the endmember induction and non-negative sparse regression using fractional semi-norm for the abundances retrieval).

In both cases, the quality of the unmixing for pixel spectrum x is evaluated by its reconstruction error (equation ( 3)).

Results

Figure 2 presents the obtained results for the semantic classes Asphalt, Vegetation and Metallic roofs. The first row of figure 2 displays the clusters obtained by the proposed strategy, where each blue spectrum depicts a local endmember obtained by LSU, the red spectrum is the cluster centroid, and the black spectrum is the corresponding endmember induced using the classical global approach. The second and third rows of figure 2 exhibit the fractional abundance maps for the classical global approach, and the proposed approach, respectively (with scales ranging from 0 (blue) to 1 (red). Finally, the bottom row displays the global scaling factors obtained by the proposed approach. Their values range from 0.5 (in blue) to 1.5 (in red). As remarked in [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF], scaling factors are only relevant if the associated fractional abundances are high enough (greater than 0.3 in practice). In the opposite situation, the contribution of the associated endmembers to the pixel spectra cannot be considered significant enough to reliably estimate their variability. In such case, scaling factors have been rounded to 1 for visualization purposes.

As it can be seen on the top row, the obtained clusters are spectrally coherent in the sense that all local endmembers grouped in the same cluster differ only from a scaling factor, which empirically validates the ELMM base assumption. Cluster centroids and classical global endmembers are also similar, up to some scaling factor, which confirms that the former can indeed be considered as global endmember instances. Nevertheless, while the obtained fractional abundance maps appear comparable for the Asphalt and Vegetation classes, which are well present across the image, it is different for the scarce Metallic roofs class. Several metallic roof endmembers have been extracted thanks to the LSU approach and clustered together, allowing to retrieve a clean abundance map, while the global approach leads to an abundance map where other structures are visible. Finally, the obtained scaling factor maps appear visually consistent, as the observed variations can be linked to the different shades within the parking lots for the Asphalt class for instance, or to the topography of the scene, as it is the case for Vegetation class (where it is possible to distinguish between trees with low scaling factors, and grass with higher scaling factors). It is even clearer when looking at the scaling factors associated to the Metallic roof class, as shown in figure 3. While the abundances of the classical global approach show some variability due to the topography of the two red roofs, this variability is clearly supported in our approach by the scaling factors, while the abundances remain relatively pure. Finally, figure 4 displays the reconstruction error maps of the proposed (figure 4a) and classical global (figure 4b) approaches. As demonstrated by figure 4c, the proposed approach globally yields lower reconstruction errors (all white pixels in figure 4c) than the classical approach. It confirms that processing LSU results within the ELMM framework allows to take advantage of both the local validity of the LMM and the global variability of endmembers to better model the hyperspectral data set.

CONCLUSION

We presented here a novel methodology to interpret LSU results within the framework defined by the ELMM. By clustering all local endmembers provided by the LSU, and defining all cluster centroids to be the ELMM reference endmembers, we showed that is was not only possible to retrieve their associated global fractional abundances, but also their scaling factors advocated by the ELMM and accounting for the endmember spectral variability (which can be due to topographic or illumination effects). The proposed method was applied to a real hyperspectral data set, and the obtained results confirmed the soundness of combining the LSU to the ELMM framework. Future work will be devoted to a exhaustive analysis of the proposed methodology over synthetic data sets, where is it possible to manually control the scaling factors and thus quantitatively assess the performance of the method. Sensitivity to the crucial choice of a clustering method will also be investigated.

Fig. 1 :

 1 Fig. 1: (a) RGB composition of the Houston hyperspectral data set, and (b) resulting segmentation composed of 396 regions.

Fig. 2 :

 2 Fig. 2: First row: obtained clusters (in blue) along with their centroids (in red) and the classical global endmembers (in black). Second and third rows: fractional abundance maps associated to the global endmembers and the cluster centroids, respectively. Bottom row: scaling factor maps obtained by the proposed methodology.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Top row: crop of the RGB image, and global abundances associated to the Metallic roof class. Bottom row: obtained scaling factors and abundances for the Metallic roof class.


G Tochon

D Pauwels

M Dalla Mura

J Chanussot

UNMIXING-BASED GAS PLUME TRACKING IN LWIR HYPERSPECTRAL VIDEO SEQUENCES

Keywords: Spectral unmixing, object tracking, hyperspectral video sequence, chemical gas plume tracking

It is now possible to collect hyperspectral video sequences (HVS) at a near real-time frame rate. The wealth of spectral, spatial and temporal information of those sequences is particularly appealing for chemical gas plume tracking. Existing state-of-the-art methods for such applications however produce only a binary information regarding the position and shape of the gas plume in the HVS. Here, we introduce a novel method relying on spectral unmixing considerations to perform chemical gas plume tracking, which provides information related to the gas plume concentration in addition to its spatial localization. The proposed approach is validated and compared with three state-of-the-art methods on a real HVS.

INTRODUCTION

Thanks to the recent advances in sensor design, it is now possible to acquire hyperspectral video sequences (HVS) at near real-time frame rates. The temporal information additionally provided by such sequences, with respect to traditional hyperspectral imagery, is particularly appealing to monitor both spectrally and spatially transient phenomena. There is however a need to adapt classical hyperspectral processing tools to handle equally the spectral, spatial and temporal information contained in HVS. Chemical gas plume tracking is a particular application that would benefit from the development of such adapted tools. As a matter of fact, this application is of interest for several real life scenarios, such as the detection of gas leaks in pipelines in order to minimize their impact on the environment [START_REF] Hinnrichs | Imaging spectrometer for fugitive gas leak detection[END_REF], or to prevent the use of chemical weapons in a context of defense and security [START_REF] Farley | Chemical agent detection and identification with a hyperspectral imaging infrared sensor[END_REF]. Most gases are in practice invisible to the human eye and respond only in a restrained portion of the long-wave infrared (LWIR) domain, hence the necessity of a fine sampling of the electromagnetic spectrum. The gas plume to monitor may in addition be very volatile and diffuse in a short lapse of time. Real-time HVS therefore appears as the most suited type of imagery for such tracking applications. State-of-the-art methods devoted to the tracking of gas plumes in HVS can be divided into three categories. Statistical methods consider the gas plume as anomalous with respect to the background and make use of anomaly detectors, see [START_REF] Manolakis | Long-Wave Infrared Hyperspectral Remote Sensing of Chemical Clouds: A focus on signal processing approaches[END_REF][START_REF] Broadwater | A Primer for Chemical Plume Detection Using LWIR Sensors[END_REF][START_REF] O'donnell | Identification and detection of gaseous effluents from hyperspectral imagery using invariant algorithms[END_REF] and references therein. Clustering-based methods aim at separating the plume from the background thanks to their spectral dissimilarities [START_REF] Merkurjev | Graph MBO method for multiclass segmentation of hyperspectral stand-off detection video[END_REF][START_REF] Qin | Detecting plumes in LWIR using robust nonnegative matrix factorization with graph-based initialization[END_REF][START_REF] Hu | Multi-class graph mumford-shah model for plume detection using the mbo scheme[END_REF]. Object tracking approaches combine temporal information and object detection methods to iteratively track the plume through the sequence [START_REF] Tochon | Gas plume detection and tracking in hyperspectral video sequences using binary partition trees[END_REF][START_REF] Tochon | Hierarchical representation of hyperspectral video sequences: application to chemical gas plume tracking[END_REF]. To the best of our knowledge however, all existing methods retrieve the position and shape of the plume in all frames of the HVS in a binary fashion, and none of them integrate the information related to its concentration. We propose to bridge this gap here by presenting a novel method for gas plume tracking that relies on spectral unmixing considerations, therefore providing this additional knowledge about the plume concentration. The rest of this paper is organized as follows: section 2 recalls the basics of spectral unmixing and introduces the notations used throughout. Section 3 features the proposed unmixingbased tracking method. Results are presented in section 4 while conclusions and perspectives are drawn in section 5.

DEFINITIONS AND NOTATIONS

Spectral unmixing

Let X ∈ R N λ ×Np be a hyperspectral data cube composed of N p pixels and N λ wavelengths. The goal of spectral unmixing is to factorize X as the product of an endmember matrix E ∈ R N λ ×d with a fractional abundances matrix Φ ∈ R d×Np , plus a noise matrix N ∈ R N λ ×Np , with d being the number of endmembers:

X = EΦ + N . (1)
Falling in the scope of the linear mixing model [START_REF] Keshava | Spectral unmixing[END_REF], each pixel spectrum x is therefore written as a linear combination of the various endmember spectra e i (corresponding to pure macroscopic elements) weighted by their corresponding fractional abundances φ i plus some noise:

x = d i=1 φ i e i + η , (2)
with the fractional abundances subjects to positivity (φ i ≥ 0 ∀i) and sum-to-one ( d i=1 φ i = 1) constraints (see [START_REF] Bioucas-Dias | Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-Based Approaches[END_REF] for an overview of spectral unmixing methods). The quality of the unmixing process for a given pixel is frequently evaluated by the root mean square error (RMSE) between the true pixel spectrum x and x = d i=1 φ i e i , reconstructed using the induced endmembers and abundances:

(x, x) = 1 N λ x -x 2 (3)

Hyperspectral video sequences

Manipulating HSV implies to handle the temporal dimension in addition to the spectral and spatial ones. For this purpose, X t will refer to the hyperspectral image acquired at time t, where the temporal axis has been discretized following the HSV frame-rate: X t+1 will therefore correspond to the frame following X t . E t = e t 1 , . . . , e t dt will denote the endmember matrix at time t (note that the number of endmembers d t may vary along the sequence). The same extension will be considered for the fractional abundances matrix Φ t .

PROPOSED UNMIXING-BASED TRACKING METHOD

Pre-processing

In a context of passive surveillance, it is realistic to assume that the recorded HVS contains plenty of frames featuring only the still background X ref (up to some minor variations due to sensor noise). The first step of the proposed unmixingbased tracking method is to extract a set of d ref reference endmembers E ref and fractional abundance Φ ref from X ref .

To overcome the noise influence, we propose in practice to perform this endmember induction and abundances retrieval over a time-averaged version of the still background.

Detection step

Considering again a passive surveillance scenario, the time where the plume appears in the sequence is likely to be a priori unknown. A detection step must therefore be performed prior to the tracking. An input frame X t either corresponds to the still background or to the plume appearing over the background in some parts of the frame. In the former case, X t should have a low reconstruction error with respect to E ref and Φ ref , while this reconstruction error should be higher in the area where the plume has appeared in the latter case. Therefore, a change detection procedure is conducted over all pixels composing the reconstruction maps (X

H 0 : ∆ t-1,t (x, x ref ) ∼ N 0, σ 2 ref H 1 : ∆ t-1,t (x, x ref ) ∼ N µ = 0, σ 2 ref (4)
where the null hypothesis H 0 corresponds to the absence of change in x between t -1 and t and the alternative hypothesis H 1 is declared when significant change has occurred. The reference variance σ 2 ref can be reliably estimated using the differences between the reconstruction maps of several background frames. The hypothesis test (4) is solved using the following Z-test statistic [START_REF] Scharf | Statistical signal processing[END_REF]:

z(x) = ∆ t-1,t (V(x), V(x ref )) σ ref / |V(x)| H1 ≷ H0 γ z ( 5
)
where ∆ t-1,t (V(x), V(x ref )) denotes the average difference reconstruction error of pixels within the neighborhood V(x) centered at x (set in practice to be a 3 × 3 neighborhood). It is known that z(x) ∼ N (0, 1) under H 0 , which allows to set the threshold γ z to achieve a predefined probability of false alarm p F A . Performing test ( 5) over all pixels of ∆ t-1,t (x, x ref ) finally leads to the binary change map C t-1,t . The plume is declared to have appeared in the sequence if C t-1,t is not empty after removal of single isolated true pixels.

Tracking step

The tracking step is launched once the gas plume has been detected in the sequence, that is when the binary change mask C t-1,t output by the detection step contains pixels where some change has been flagged. Likewise classical object tracking methods (see [START_REF] Yilmaz | Object tracking: A survey[END_REF] for a review), the proposed approach is decomposed in two inner steps addressed sequentially.

Motion prediction step

The goal of the motion prediction step is to estimate an area in the current frame X t where the gas plume can be found with a high probability. Here, the change detection procedure previously described in section 3.2 is applied to X t-1 and X t , to produce a binary change map C t-1,t featuring areas that have changed between the two consecutive frames. More particularly, C t-1,t describes the regions that have either been left by the plume, or invaded by it between t -1 and t. Denoting by O t-1 the binary mask of the plume position in X t-1 , it is possible to obtain an estimate position of the plume in X t by

Ôt = O t-1 C t-1,t (6)
where stands for the binary XOR operation, following [START_REF] Tochon | Gas plume detection and tracking in hyperspectral video sequences using binary partition trees[END_REF]. Note that Ôt = C t-1,t for the frame where the plume appears, as O t-1 is empty in this case. The proposed motion prediction step is illustrated by the workflow presented by figure 1.

Xt Xt-1 (Xt, Xref ) (Xt-1, Xref ) ∆ t-1,t (X, Xref ) Ct-1,t ⊕ Ôt Ot-1
Fig. 1: Workflow of the proposed motion prediction step.

Matching step

The goal of the matching step is to retrieve the position of the plume O t in the current frame X t given the estimation Ôt provided by the motion prediction step. Here, a set of d ref + 1 endmembers E t is generated using some endmember induction algorithm (EIA), where d ref is the number of reference endmembers used to describe the static background.

Ideally, the EIA should generate one endmember associated to the gas plume signature, and the others corresponding to the background. Being in an unsupervised context where the theoretical gas plume signature is unknown, the corresponding endmember e t p is identified as the one that falls in the estimated plume area Ôt . If several or no endmembers are induced within Ôt , then the EIA can be run again until exactly one endmembers falls in Ôt , or e t p can be enforced as the spectrum of the pixel having the highest reconstruction error within Ôt in (X t , X ref ). Once the current endmembers E t have been generated, their corresponding fractional abundances Φ t are retrieved. Following, a classification map Θ t of X t is issued using a maximum voting rule on the abundances Φ t . More specifically, let φ t = φ t 1 , . . . , φ

assigns each pixel x t to the class corresponding to the endmember whose fractional abundance is the highest in x t . Following, the classification map Θ t is used to produce the binary map O t featuring only the pixels classified in the gas plume class. This binary map is applied over the abundance map associated to the gas plume endmember e t p in order to finally obtain the data of the position, shape, and concentration of the plume O t in X t . This whole matching process is summarized by the workflow in figure 2.

EXPERIMENTS

Data set

The data set used to validated the proposed unmixing-based tracking methodology described in the previous section 3 was provided by John Hopkins Applied Physics Laboratory. The video sequence was acquired in 2006 at the Dugway Proving Ground in Utah (USA) with a FIRST sensor [START_REF] Farley | Performance of the first: a longwave infrared hyperspectral imaging sensor[END_REF] at a frame-rate of 0.2 Hz. Each frame is a data cube composed of 128 × 320 pixels in the spatial domain, and 129 spectral channels spanning 7.81 µm to 11.97 µm in wavelength. The HVS features the explosive release of an acetic acid canister and the resulting sequence is composed of 30 frames, with the gas release occurring at the 11 th frame. Ground truth data was created in order to quantitatively evaluate the performances of the proposed method. A false color composition of each frame was first derived, and two classes (corresponding to strongly concentrated and more diffused areas within the plume) were then carefully delineated. A false color composition of the 20 th frame along with its ground truth map are displayed by figure 3.

Experimental methodology

The critical parameter to operate the proposed unmixing-based tracking methodology is the definition of the EIA, as inconsistent endmembers would not yield proper results. Among the broad literature of EIA [START_REF] Bioucas-Dias | Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-Based Approaches[END_REF], we choose to work with the VCA algorithm [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] as it already proved to efficiently unmix LWIR hyperspectral images with faint gas signatures [START_REF] Henrot | Dynamical spectral unmixing of multitemporal hyperspectral images[END_REF]. Fractional abundance are then estimated using the SUNSAL [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF] method. Following [START_REF] Broadwater | A Primer for Chemical Plume Detection Using LWIR Sensors[END_REF], the number of reference endmembers d ref to described the background is set to 4. It is assumed that the first 5 frames of the sequence are known not to contain the plume. Consequently, the reference variance σ 2 ref is estimated over the first 4 differences of reconstruction error maps. The Z-test statistic ( 5) is operated under a probability of false alarm set to p F A = 0.01. Due to the stochasticity of the VCA, the proposed method is run 20 times, and the average percentages of correct detections (amount of ground truth pixels correctly detected) and false detections (amount of ground truth pixels wrongly detected) with respect to the ground truth are computed for all frames of the sequence.

Results

Figure 5 shows the average percentage of correct and false detections for the proposed unmixing-based method (with standard deviation), as well as for three state-of-the-art approaches, namely the AMSD [START_REF] Broadwater | A Primer for Chemical Plume Detection Using LWIR Sensors[END_REF] (statistical-based), RNMF [START_REF] Merkurjev | Graph MBO method for multiclass segmentation of hyperspectral stand-off detection video[END_REF] (clustering-based) and BPT [START_REF] Tochon | Hierarchical representation of hyperspectral video sequences: application to chemical gas plume tracking[END_REF] (object detection based) methods. As it can be seen in figure 5a, the proposed approach performs equally well than the RMNF approach, and outperforms both AMSD and BPT-based methods for most frames of the sequence. While the release instant of the plume is correctly identified, the proposed method fails at properly detecting the gas plume in the early frames #11 and #12. The reason is that the EIA is not able to automatically induce an endmember in the estimated area Ôt , which is then forced to be the pixel with the highest reconstruction error. Due to the explosive release of the canister containing the gas, the enforced plume endmember is spectrally aberrant and unable to correctly describe the gas plume, hence the failing of the subsequent classification step. Nonetheless, the average percentage of correct detections for the other frames of the sequence remains around 80% while the AMSD and BPT methods suffer a drop as the diffusion is becoming more severe. It is also worth adding that the RMNF method processes all frames of the sequence in a stacked manner, therefore being incompatible with a real-time implementation scenario where incoming frames would have to be processed on the fly. Figure 5b shows the percentage of false detections, corresponding to the pixels whose spectra are supposedly not corrupted by the plume but which are nonetheless detected by the various methods. The y-axis has been cut for a better visualization. Apart for the AMSD approach that produces a significant amount of false detections (notably before the release of the plume), all three other competing approaches correctly identify the release instant, and their number of false detections remain relatively low across the sequence. The previous observations are confirmed by the qualitative analysis of figure 4. While the proposed approach is not able to properly detect the plume for the early frame #12, the retrieved shape and concentration appear visually very consistent for all the other frames of the sequence.

CONCLUSION

We presented here a novel method to perform the detection and tracking of a diffusing gas plume in HVS. Unlike classical state-of-the-art methods, the proposed approach is not only able to retrieve the position and shape of the plume in all frames of the sequence, but also provides information related to its concentration thanks to the integration of a spectral unmixing process within the tracking procedure. The proposed method was investigated over a real HVS and compared with three state-of-the-art approaches, and showed to perform equally or better while providing a richer information. Extending the proposed approach to accommodate non linear spectral mixing models is part of our research perspectives, as well as its optimized implementation to make it operable under real-time processing constraints.

Fig. 2 :

 2 Fig. 2: Workflow of the proposed matching step.

Fig. 3 :

 3 Fig. 3: (left) False color composition of the 20 th frame of the aa12 HVS and (right) corresponding ground truth data.

Fig. 4 :

 4 Fig. 4: Top row: frames #12, 16, 20, 24 and 28 of the sequence. Bottom row: corresponding detection O t .

Fig. 5 :

 5 Fig. 5: Percentage of (a) correct and (b) false detections with respect to the ground truth.

  the plume has contaminated x t , since E ref and Φ ref alone cannot properly reconstruct the plume signature. We can therefore design the following hypothesis test:

t-1 , X ref ) and (X t , X ref ). More specifically, let (x t-1 , x ref ) and (x t , x ref ) be the reconstruction errors of a given pixel spectrum at t -1 and t, and let ∆ t-1,t (x, x ref ) = (x t , x ref ) -(x t-1 , x ref ) be the difference of such errors. If no change has occurred between t -1 and t within pixel x, this spectrum should have relatively close reconstruction error values, while their difference should be higher if

ACKNOWLEDGMENTS

We thank Andrea Bertozzi, Department of Mathematics, University of California, Los Angeles (UCLA), for her collaboration on the dataset provided by the US Defense Threat Reduction Agency and the National Science Foundation through NSF grant DMS-1118971. This work was also partially funded through DGA, under grant 2015-60-001200.470.75.01.


Alexis Hannart
email: alexis.hannart@cima.fcen.uba.ar

Alberto Carrassi

Marc Bocquet

Michael Ghil

Philippe Naveau

Manuel Pulido

Juan Ruiz

Pierre Tandeo

A Hannart Ifaeci

DADA: Data Assimilation for the Detection and Attribution of weather-and climate-related events

similation procedures that are run on a daily basis to update the meteorological model with new atmospheric observations; hence, the proposed methodology can take advantage of the powerful computational and observational capacity of weather forecasting centers. We explain the theoretical rationale of this approach and sketch the most prominent features of a "data assimilation based detection and attribution" (DADA) procedure. The proposal is illustrated in the context of the classical three-variable Lorenz model with additional forcing. Several theoretical and practical research questions that need to be addressed to make the proposal readily operational within weather forecasting centers are finally laid out.

Keywords Event attribution • Data assimilation • Causality theory • Modified Lorenz model 1 Background and motivation A significant and growing part of climate research studies the causal links between climate forcings and observed responses. This part has been consolidated into a research topic known as detection and attribution (D&A). The D&A community has increasingly been faced with the challenge of generating causal information about episodes of extreme weather or unusual climate conditions. This challenge arises from the needs for public dissemination, litigation in a legal context, adaptation to climate change or simply improvement of the science associated with these events [START_REF] Stott | Attribution of weather and climate-related events[END_REF].

The approach widely used so far to in D&A was introduced one decade ago by M.R. Allen and colleagues [START_REF] Allen | Liability for climate change[END_REF][START_REF] Stone | The end-to-end attribution problem: from emissions to impacts[END_REF] and it originates from best practices in epidemiology [START_REF] Greenland | Measures of effect and measures of association[END_REF]. In this approach, one evaluates the extent to which a given external climate forcingsuch as solar irradiation, greenhouse gas (GHG) emissions, ozone or aerosol concentrations -has changed the probability of occurrence of an event of interest.

For this purpose, one compares the probability of occurrence of said event in an ensemble of model simulations representing the observed climatic conditions, which simulates the actual occurrence probability in the real world, with the occurrence probability of the same event in a parallel ensemble of model simulations, which represent an alternative world. The former world is referred to as factual, the latter as counterfactual: it is the latter that might have occurred had the forcing of interest been absent.

Denoting by p 1 and p 0 the probabilities of the event occurring in the factual world and in the counterfactual world respectively, the so-called fraction of attributable risk (FAR) is then defined as FAR= 1-p 0 /p 1 . The FAR has long been interpreted as the fraction of the likelihood of an event which is attributable to the external forcing. Over the past decade, most causal claims have been following from the FAR and its uncertainty, resulting in statements such as "It is very likely that over half the risk of European summer temperature anomalies exceeding a threshold of 1.6 • C is attributable to human influence." [START_REF] Stott | Human contribution to the European heatwave of 2003[END_REF]. [START_REF] Hannart | Counterfactual causality theory for the attribution of weather and climate-related events[END_REF] have recently shown that, under realistic assumptions, the FAR may also be interpreted as the so-called probability of necessary causation (PN) associated -in a complete and self-consistent theory of causality [START_REF] Pearl | Causality: Models, Reasoning and Inference[END_REF] -with the causal link between the forcing and the event. The FAR thus corresponds to only one of the two facets of causality in such a theory, while the probability of sufficient causation (PS) is its second facet. In this setting,

PN = 1 - p 0 p 1 , (1a)
PS = 1 - 1 -p 1 1 -p 0 , ( 1b
) PNS = p 1 -p 0 , ( 1c
)
where PNS is the probability of necessary and sufficient causation. [START_REF] Pearl | Causality: Models, Reasoning and Inference[END_REF] provides rigorous definitions of these three concepts, as well as a detailed discussion of their meanings and implications. It can be seen from Eqs.

(1) that causal attribution requires to evaluate the two probabilities, p 0 and p 1 , and not just one of them. Doing so is, therefore, the central methodological question of D&A for weather and climate-related events. So far, most case studies have used large ensembles of climate model simulations in order to estimate p 1 and p 0 based on a variety of methods, in particular based on statistical extreme value theory (EVT). However, this general approach has a very high computational cost and is difficult to implement in a timely and systematic way. As recognized by [START_REF] Stott | Attribution of weather and climate-related events[END_REF], this remains an open problem: "the overarching challenge for the community is to move beyond research-mode case studies and to develop systems that can deliver regular, reliable and timely assessments in the aftermath of notable weather and climate-related events, typically in the weeks or months following (and not many years later as is the case with some researchmode studies)". For instance, the weather@home system [START_REF] Massey | weather@home -development and validation of a very large ensemble modelling system for probabilistic event attribution[END_REF], or the system proposed by [START_REF] Christidis | A New HadGEM3-A-Based System for Attribution of Weather-and Climate-Related Extreme Events[END_REF], aim at meeting those requirements within the conventional ensemble-based approach. Ongoing research aiming towards the development of such a system also include the CASCADE project (Calibrated and Systematic Characterization, Attribution and Detection of Extremes, U.S. Department of Energy, Regional and Global Climate Modeling program).

The purpose of this article is to introduce a new methodological approach that addresses the latter overarching operational challenge. Our proposal relies on a class of powerful statistical methods for interfacing high-dimensional models with large observational datasets. This class of methods originates from the field of weather forecasting and is referred to as data assimilation (DA) (Bengtsson et al., 1981;[START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF][START_REF] Talagrand | Assimilation of observations, an introduction[END_REF].

Section 2 explains the rationale of the approach proposed herein, presents a brief overview of DA, and outlines the most prominent technical features of a "data assimilation-based detection and attribution" (DADA) approach. Section 3 illustrates the proposal by implementing it on a version of the classical Lorenz convection model (Lorenz, 1963, L63 hereafter) subject to an additional constant force. Finally, in Section 4, we discuss the main strengths and limitations of the DADA approach, and highlight several theoretical and practical research questions that need to be addressed to make it potentially operational within weather forecasting centers in a near future.

Method description

General rationale

The rationale for addressing causal attribution of climaterelated events based on DA concepts and methods can be outlined in three steps. To do so briefly and clearly, we need to introduce some notation.

Let y t denote the d-dimensional vector of observations at discrete times {t = 0, 1, . . . , T }. Here, y = {y t : 0 ≤ t ≤ T } corresponds, for instance, to the full set of all available meteorological observations over a time interval covering the event of interest, no matter the diversity and source of the data; typically, the latter include ground station networks, satellite measurements, ship data, and so on, cf. (Bengtsson et al., 1981, Preface, Fig. 1) or [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF], Fig. 1). In the present probabilistic D&A context, the observed trajectory y is viewed as a realization of a random variable denoted Y = {Y t : 0 ≤ t ≤ T }, i.e. there exists an ω ∈ Ω such that Y(ω) = y -where Ω denotes the sample space of all possible outcomes and encompasses observational error, as well as internal variability.

In event attribution studies, it is recognized that defining the occurrence of an event, i.e. selecting a subset F ⊂ Ω, depends on a rather arbitrary choice. Yet this choice has been shown to greatly affect causal conclusions [START_REF] Hannart | Counterfactual causality theory for the attribution of weather and climate-related events[END_REF]. For instance, a generic and fairly loose event definition is arguably prone to yield a low threshold of evidence with respect to both necessary and sufficient causality while, on the other hand, a tighter and more specific event definition is prone to yield a stringent threshold for necessary causality but a reduced one for sufficient causality.

Indeed, it is quite intuitive that many different factors should usually be necessary to trigger the occurrence of a highly specific event and conversely, that no single factor will ever hold as a sufficient explanation thereof. For the class of unusual events at stake in D&A, where both p 0 and p 1 are very small, we arguably lean towards specific definitions that inherently result in few sufficient causal factors or none. This conclusion immediately follows from Eq. (1b), which yields PS 0 when both p 0 and p 1 are very small.

Usually, an event occurrence is defined in D&A based on an ad hoc scalar index φ(Y) exceeding a threshold u, i.e. p i = P (φ(Y) ≥ u); from now on, we associate i = 0 with the counterfactual and i = 1 with the factual world. While this definition may be already quite restrictive for u large, it is a defensible strategy to restrict the event definition even further: this may slightly reduce an already negligible PS but in return may potentially increase PN by a greater amount; one thus expects to gain more than one loses in this trade-off. In particular, this will be the case if additional features, not accounted for in φ(Y), can be identified that will allow one to further discriminate between the two worlds.

In any case, a central element of our proposal is to follow this strategy in its simplest possible form, by using the tightest occurrence definition i.e. the singleton {ω ∈ Ω | Y(ω) = y}. Note that the latter singleton has probability zero in both worlds because the probability density function (PDF) f (Y(ω)) of Y can be assumed, in general, to be continuous, i.e. to contain no singular δ-functions.

Consider, however, the paradox that arises from taking the limit h → 0 for the set {ω ∈ Ω | Y(ω) -y ≤ h}. This set has non-zero probability for h arbitrarily small but positive while, in the limit,

PN = 1 - f 0 (y) f 1 (y) , PS = 0, (2)
where f i denotes the PDF of Y in world i. Equation ( 2) thus shows that, while the probabilities of occurrence of our singleton event in both worlds are null, its associated probability of necessary causation is still positive -but its probability of sufficient causation is always zero. Our proposal thus intentionally sacrifices evidence of sufficiency, in the hope of maximizing the evidence of necessity.

Our betting on the singleton set is thus justifiable already based on the above theoretical considerations. This choice, moreover, is motivated by having a highly simplifying implication from a practical standpoint. Evaluating the PDF of Y at a single point Y = y is indeed, under many circumstances, considerably easier than evaluating the probability P (φ(Y) ≥ u) required in the conventional approach.

To illustrate this point, let Y be for instance a dvariate autoregressive process defined by Y t+1 = AY t + w t , where w t is an i.i.d. noise having known PDF g(•) and where A has the usual properties that insure stationarity [START_REF] Gardiner | Handbook of Stochastic Methods for Physics[END_REF]. We then have:

f (y) = T t=1 g(y t -Ay t-1 )π(y 0 ) , (3a)
P (φ(Y) ≥ u) = φ(y)≥u T t=1 g(y t -Ay t-1 ) × π(y 0 )dy 1,0 . . . dy d,0 . . . dy d,T , (3b)
with π(•) the prior PDF on the initial state Y 0 . Equation (3a) shows that f (y) can be easily computed using a closed-form expression, while P (φ(Y) ≥ u) in Eq. ( 3b) is an integral on d × T + 1 dimensions which must instead be evaluated by using, for instance, a computationally quite costly Monte-Carlo (MC) simulation. Figure 1 illustrates this situation by showing the details of the latter MC evaluation for a scalar AR [START_REF] Kondrashov | Log-normal Kalman filter for assimilating phase-space density data in the radiation belts[END_REF] process (panel a, when based on a standard EVT application, as well as its associated accuracy (panels b and c), and the computational cost as the MC sample size n varies (panel d); the latter cost is much larger than the one of applying the DADA approach. This simple example confirms the large computational discrepancy between the two approaches.

The reason for the discrepancy is quite simple: evaluating the conventional probability requires integrating a PDF over a predefined domain, instead of a one-off evaluation at a single point. Because both the domain of integration and the PDF may have potentially complex shapes, one cannot expect, in general, that the requisite integral be amenable to analytical treatment. Hence numerical integration is the default option: no matter how efficient an integration scheme one applies, it will require evaluating the PDF at many points and is thus as many times more costly computationally than just evaluating f (y) at a single point.

This being said, it is not always straightforward to obtain the PDF of Y. This is the case, for instance, for the wide class of statistical models referred to as Hidden Markov Models (HMMs); in fact, HMMs [e.g., (Ihler et al., 2007, and references therein)] are often relevant in the present context to describe Y.

More precisely, assume that the event of interest can be represented by a large numerical model which N -dimensional state vector at time t is denoted X t . The dynamics of the state vector is given by:

X t+1 = M(X t , F t ) + v t , ( 4
)
where M is the model operator, v t is a stochastic term representing modeling error, and F t is a known, prescribed forcing that is external to the model. In the present context, it is precisely the forcing term F = (F t ) T t=0 that is under causal scrutiny. Further, assume that our observations Y t can be mapped to the state vector X t at any time t, i.e.

Y t = H(X t ) + w t (5)
where H is the so-called observation or forward operator and w t is a stochastic term representing observational error.

Denoting by F (i) the value of the forcing in the world i, using the shorthand M i (x t ) = M(x t , F (i) t ) and denoting by M i the HMM associated with H and M i , the problem of interest here is thus to derive:

f 0 (y) = f (y | M 0 ) and f 1 (y) = f (y | M 1 ) , (6)
where f 0 (y) and f 1 (y) should be interpreted as the likelihoods of the observation y in the counterfactual and factual models, respectively. Finally getting to our point, one can view DA methods as a class of inference methods designed for the above HMM setting. Actually, [START_REF] Ihler | Graphical models for statistical inference and data assimilation[END_REF] already formulated both DA and HMMs within the broader class of graphical models for statistical inference.

While inferring the unknown state vector trajectory X, given the observed trajectory y, is clearly the main focus of DA, the likelihood f (y) can also be obtained as a side product thereof, as we will immediately clarify below. Therefore, with DA able to derive the two likelihoods f 0 (y) and f 1 (y), and the latter two being the keys to causal attribution in our approach, one should be capable of moving towards near-real-time, systematic causal attribution of weather-and climate-related events.

Brief overview of data assimilation

DA was initially developed in the context of numerical weather forecasting, in order to initialize the model's state variables X based on observations y that are incomplete, diverse in nature, unevenly distributed in space and time, do not necessarily match the model's state variables, and are contaminated by measurement error (Bengtsson et al., 1981;[START_REF] Talagrand | Assimilation of observations, an introduction[END_REF]. Over the past decades, those methods have grown out of their original application field to reach a wide variety of topics in geophysics such as oceanography [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF], atmospheric chemistry, geomagnetism, hydrology, and space physics, among many other areas [START_REF] Robert | Comparison of reduced-order sequential, variational and hybrid data assimilation methods in the context of a Tropical Pacific ocean model[END_REF]Cosme et al., 2010;[START_REF] Kondrashov | Log-normal Kalman filter for assimilating phase-space density data in the radiation belts[END_REF][START_REF] Bocquet | Parameter-field estimation for atmospheric dispersion: application to the Chernobyl accident using 4D-Var[END_REF][START_REF] Martin | Status and future of data assimilation in operational oceanography[END_REF].

DA is already playing an increasing role in the climate sciences, having being applied, for instance, to initialize a climate model for seasonal or decadal prediction [START_REF] Balmaseda | Ocean initialization for seasonal forecasts[END_REF], to constrain a climate model's parameters [START_REF] Kondrashov | Data assimilation for a coupled ocean-atmosphere model. Part II: Parameter estimation[END_REF][START_REF] Ruiz | Estimating model parameters with ensemble-based data assimilation: A review[END_REF], to infer carbon cycle fluxes from atmospheric concentrations [START_REF] Chevallier | On the parallelization of atmospheric inversions of CO2 surface fluxes within a variational framework[END_REF], or to reconstruct paleoclimatic fields out of sparse and indirect observations [START_REF] Bhend | An ensemble-based approach to climate reconstructions Clim[END_REF][START_REF] Roques | Parameter estimation for energy balance models with memory[END_REF]. In the context of D&A, [START_REF] Lee | Evaluation of proxy-based millennial reconstruction methods[END_REF] actually tested a DA-like approach to include the effects of the various forcings over the last millennium, in addition to other paleoclimate proxy data, in combined climate reconstruction and detection analysis. The present work thus follows and further strengthens a general trend in climate studies.

Methodologically speaking, DA methods are traditionally grouped into two categories: sequential and variational (Ide et al., 1997, and references therein). In the sequential approach [START_REF] Ghil | Applications of estimation theory to numerical weather prediction[END_REF], the state estimate and a suitable estimate of the associated error covariance matrix are propagated in time until new observations become available and are used to update the state estimate. In practice, the evolution of the system of interest is retrieved -like in earlier, typically much smaller-dimensional applications [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Jazwinski | Stochastic and Filtering Theory[END_REF]Gelb, 1974) -through a sequence of prediction and analysis steps. In the variational approach, on the other hand, one seeks the system trajectory that best fits all the observations distributed within a given time interval (Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects[END_REF][START_REF] Ide | Unified notation for data assimilation: Operational, sequential and variational[END_REF][START_REF] Bocquet | Parameter-field estimation for atmospheric dispersion: application to the Chernobyl accident using 4D-Var[END_REF]. Here, we concentrate on the sequential approach, but the two approaches are complementary and the choice of method depends on the specifics of the problem at hand [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF][START_REF] Ide | Unified notation for data assimilation: Operational, sequential and variational[END_REF][START_REF] Talagrand | Assimilation of observations, an introduction[END_REF].

Abundant literature is available on DA and on Kalmantype filters. [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] first presented the solution in discrete time for the case in which both the dynamic evolution operator M in Eq. 4 and the observation operator H in Eq. 5 are linear, and the errors are Gaussian. Under these assumptions, the state-estimation problem for the system given by Eqs. (5, 4) has an exact solution given by the following sequential Kalman filter (KF) equations:

x a t = x f t + K(y t -Hx f t ) , (7a)
P a t = (I -KH)P f t , (7b)
x f t+1 = Mx a t , (7c)
P f t+1 = MP a t M + Q . ( 7d
)
where denotes the transpose operation. Here Eqs. (7a) and (7b) are referred to as the analysis step and denoted by a superscript a, while the forecast step is given by Eqs. (7c) and (7d), and is denoted by a superscript f [START_REF] Ide | Unified notation for data assimilation: Operational, sequential and variational[END_REF]. The vector x a t and the matrix P a t are the mean and covariance of

X t conditional on (Y 1 , ..., Y t ) = (y 1 , ..., y t ); K = P f t H (HP f t H + R) -1
is the so-called Kalman gain matrix; while Q and R are the covariances associated with v t and w t , respectively. Following Wiener (1949), one distinguishes between filtering, in which x a t and P a t are conditioned only on the previous and current observations (y 0 ,..., y t ), and smoothing, in which they are conditioned on the entire sequence, 0 ≤ t ≤ T . Furthermore, the sequential algorithm needs to be initialized at time t = 0 with x f 0 and P f 0 , which thus represent the a priori mean and covariance of X 0 , respectively, and have to be prescribed by the user.

The likelihood function f (y), which is of primary importance for DADA, also has an exact expression under the above linearity and Gaussianity assumptions (Tandeo et al., 2014), given by:

f (y) = T t=0 (2π) -d 2 |Σ t | -1 2 × exp - 1 2 (y t -Hx f t ) Σ -1 t (y t -Hx f t ) , (8)
with Σ t = HP f t H + R. The proof of Eq. ( 8) is provided in the Appendix, and f (y) is typically computed by taking the logarithm of this equation to turn the product on the right-hand side into a sum.

It follows from the above that, once the observations y t have been assimilated on the interval 0 ≤ t ≤ T , the necessary ingredients x f t and P f t in Eq. 8 are available and thus calculating f (y) is both straightforward and computationally inexpensive. The fundamental connections between this calculation, the HMM context, and Bayes theorem are further clarified in the Appendix.

Many difficulties arise in applying the simple ideas outlined here to geophysical models, which are typically nonlinear, have non-Gaussian errors and are huge in size [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF]. Most of these difficulties have been addressed by improving both sequential and variational methods in several ingenious ways [START_REF] Bocquet | Beyond Gaussian statistical modeling in geophysical data assimilation[END_REF][START_REF] Kondrashov | Log-normal Kalman filter for assimilating phase-space density data in the radiation belts[END_REF].

In particular, the Ensemble Kalman Filter (EnKF; Evensen, 2003)-in which the uncertainty propagation is evaluated by using a finite-size ensemble of trajectories -is now operational in numerical weather and oceanic prediction centers worldwide; see e.g. Sakov et al. (2013); [START_REF] Houtekamer | Higher Resolution in an Operational Ensemble Kalman Filter[END_REF]. The EnKF is a convenient approximate solution to the filtering problem in a nonlinear, large-dimensional context. We simply note here that it can also be applied to obtain an approximation of the likelihood f (y) by substituting the approximate sequence {(x f t , Pf t ) : t = 0, . . . , T } that the EnKF produces into Eq. 8. This strategy is illustrated immediately below in the context of the L63 convection model subject to an additional constant force.

3 Implementation within the modified L63 model

The modified model and its two worlds

A simple modification [START_REF] Palmer | A non-linear dynamical perspective on climate prediction[END_REF] of the L63 system [START_REF] Lorenz | Deterministic non-periodic flow[END_REF] has been extensively used for the purpose of illustrating methodological developments in both DA and D&A [e.g. [START_REF] Carrassi | Model error and variational data assimilation: A deterministic formulation[END_REF][START_REF] Stone | The end-to-end attribution problem: from emissions to impacts[END_REF]]. In the nonlinear, coupled system of three ordinary differential equations (ODEs) for

x, y and z below,

dx dt = σ(y -x) + λ i cos θ i , dy dt = ρx -y -xz + λ i sin θ i , dz dt = xy -βz (9)
the time-constant forcing terms in the x-and y-equation represent, in fact, an addition to the forcing hidden in the original L63 model. The latter forcing is revealed by a well-known linear change of variables, in which x and y are left unchanged and z → z + ρ + σ [START_REF] Lorenz | Deterministic non-periodic flow[END_REF].

In the new variables, the model of Eq. ( 9) will take the canonical form of a forced-dissipative system (Ghil and Childress, 1987, Sec. 5.4), with an extra forcing term -β(ρ + σ) in the z-equation, just like the original L63 model.

Here λ i is the intensity of the additional forcing and θ i is its direction in world i = 0, 1: i.e., λ 0 = 0 represents a counterfactual world with no additional forcing, while λ 1 = 0. We take the parameters (σ, ρ, β) to equal their usual values (10, 28, 8/3) that yield the well-known chaotic behavior, and the (nondimensional) time unit t is interpreted as equaling days.

The ODE system given by ( 9) is discretized by using ∆t = 0.01 and t refers hereafter to the number of time increments ∆t. This system is then turned into one of stochastic difference equations [S∆Es: Arnold (2003); [START_REF] Chekroun | Stochastic climate dynamics: Random attractors and time-dependent invariant measures[END_REF]] by adding an error term v t assumed to be Gaussian and centered with covariance Q = σ 2 Q I, where I is the 3 × 3 identity matrix. Furthermore, we assume that all three coordinates (x, y, z) of the state vector are observed, i.e. that H = I, and that the measurement error term w t is also Gaussian and centered, with covariance R = σ 2 R I. Recalling the notation introduced in Sec 2a, we associate a label ω ∈ Ω with each realization of the pair of random processes (v t , w t ) that drive the model given by Eq. ( 9) and perturb its observations, respectively.

The S∆E system defined above is stationary, i.e. the PDF of the state vector x t depends neither on t nor on x 0 after a sufficiently long time t. This PDF can be obtained as the (numerical) solution of the Fokker-Planck equation associated with Eq. ( 9), and it is the mean over Ω of the sample measures obtained for each realization ω of the noises v t and w t (Chekroun et al., 2011, and references therein). Each sample measure is supported on a random attractor that may have very fine structure and be time-dependent (Chekroun et al., 2011, Figs. 1-3 and supplementary material), but the PDF is supported smoothly, in the counterfactual world in which λ 0 = 0, on a "thickened" version of the fairly well-known strange attractor of the original L63 model.

In the factual world in which λ 1 = 0, the nature of the PDF is quite similar, but its exact shape is af-fected by the parameters (λ 1 , θ 1 ) of the forcing. In both worlds, the PDFs can be estimated, for instance, by using kernel density estimation applied to ensembles of simulations obtained for either forcing. In Figs. 2a,b, we plot the projections of both PDFs onto the plane associated with the greatest variance in the factual PDF. The difference between the two PDFs is shown in Fig. 2c; it emphasizes the existence of an area of the state space (represented in white), which is more likely to be reached in the factual world than in the counterfactual one.

Next, we define an event to occur for the sequence {y t : t = 0, . . . , T } if the scalar product φ y t between the unit vector φ in the direction φ and y t , i.e. the projection of y t onto the direction φ, exceeds u for some 0 ≤ t ≤ T , where φ is a specified direction and u is a threshold chosen based on φ so that p 1 = 0.01. Figure 2d shows a selection of sequences from both worlds in which an event did occur, where φ was chosen to be the leading direction in the projection plane.

For this choice of φ, the trajectories associated with event occurrence happen to all lie in the area of the state space which is more likely to be reached in the factual world than in the counterfactual one. Accordingly, the probability of the event in the former is found to be higher than in the latter, i.e. p 1 > p 0 , and the occurrence of an event {max {0≤t≤T } φ y t ≥ u} is thereby informative from a causal perspective, i.e. the associated probabilities of necessary and sufficient causation are positive.

Figure 2d also shows that the trajectories associated with the event in the two worlds -counterfactual (green) and factual (red) --appear to have slightly distinct features: the red trajectories are shifted towards higher values in the second direction, of highest-but-one variance. Such distinctions might help discriminate further between the two worlds in the DADA framework.

DADA for the modified L63 model

The DADA procedure is illustrated in Fig. 3. We plot in panel (a) a trajectory of the state vector x t simulated under factual conditions, i.e. in the presence of the additional forcing (black solid line), along with the observations {y t : 0 ≤ t ≤ T } (gray dots), with T = 400. The EnKF is used to assimilate these observations into a factual model (i = 1) that thus matches the true model M = M 1 = M(λ 1 , θ 1 ) used for the simulation: a reconstructed trajectory is obtained from the corresponding analyses x a t (red solid line in panel (a)), cf. Eqs. ( 7), and the likelihoods f 1 (y t ) (red solid line in panel (c)) are obtained by application of Eq. ( 8), respectively.

Next, the assimilation is repeated in the counterfactual model (i = 0, i.e. λ = 0) to obtain a second analysis of the trajectory, from the same observations; see green solid line in panel (a), for T = 400. The corresponding likelihoods f 0 (y t ) are shown in panel (c) as a green solid line. Comparing the trajectories of the two analyses in Fig. 3a shows that, even though the counterfactual analysis (green line) uses the same data as the factual analysis (red line), the former lies closer to the true trajectory (black line).

The local discrepancies between the trajectories estimated in the two worlds appear to be rather small at first glance, cf. panel (a), and so are the instantaneous differences between the associated factors on the right-hand side of Eq. ( 8); the latter are shown as gray rectangles in panel (c) of the figure. Still, the evidence in favor of the factual world accumulates as the time t over which the two trajectories differ, albeit by a small amount, lengthens. This cumulative difference in evidence, log f 0 (y t ) -log f 1 (y t ), is reflected by a growing gap between the two curves, red and green, in panel (c), and by an associated high mean growth over time of the probability PN of necessary causation, cf. the black solid line in panel (d).

In order to evaluate more systematically its performance and robustness compared to the conventional FAR approach, the DADA procedure was applied to a large sample of sequences y t of length T = 20 simulated under diverse conditions. The sample explored all possible combinations of the triplet of parameters (λ 1 , σ Q , σ R ), with ten equidistributed values each, for a total of 10 3 combinations; the ranges were 0 ≤ λ 1 ≤ 40, 0.1 ≤ σ Q ≤ 0.5 and 0.1 ≤ σ R ≤ 1.0, respectively, with θ 1 = -140 • . For each combination of (λ 1 , σ Q , σ R ), ten directions φ were randomly generated and u was defined based on φ as in Sec. 3a above, so as to achieve p 1 ≥ 0.01.

In order to estimate the corresponding conventional probabilities p 0 and p 1 of the associated event defined as {max {0≤t≤T } φ y t ≥ u}, n = 50 000 sequences y t of length T = 20 were simulated, by using a single sequence of length nT = 10 6 and splitting it into n equal segments. Probabilities p 0 and p 1 were then directly estimated from empirical frequencies because the high value of n here did not require the use of the EVT extrapolation normally used for smaller n.

For each quintuplet of parameter values (λ 1 , σ Q , σ R ; φ, u), one hundred sequences of observations {y t : 0, . . . , T = 20} were generated with a proportion p 1 /(p 1 + p 0 ) being simulated from the factual world and a proportion p 0 /(p 1 + p 0 ) from the counterfactual one. All sequences were treated with the DADA procedure -by applying DA to the synthetic observations according to Eqs. ( 7a)-(7d) -and then Eq. ( 8) to obtain f 0 (y) and f 1 (y) from the reconstructed tra-jectories. The a priori mean and covariance x f 0 and P f 0 required as inputs to the DADA procedure were those associated with the PDF of the attractor, given the forcing conditions (λ 1 ∈ [0, 40], θ 1 = -140 • ) assumed for each assimilation experiment. As a result, two probabilities PN of necessity are finally obtained for each sequence y t , PN p = 1 -p 0 /p 1 for the conventional approach and PN f = 1 -f 0 (y)/f 1 (y) for the DADA approach.

We next wish to evaluate under various conditions how well the two probabilities PN p and PN f perform with respect to discriminating between the factual and counterfactual forcings. Consider a simple discrimination rule whereby a trajectory y t is identified as factual for PN exceeding a given threshold, and as counterfactual otherwise. The so-called receiver operating characteristic (ROC) curve plots the rate of true positives as a function of the rate of false positives obtained when varying the threshold in a binary classification scheme from 0 to 1; it thus gives an overall visual representation of the skill of our PN as a discriminative score.

The [START_REF] Gini | Measurement of inequality of incomes[END_REF] index G was originally introduced as a measure of statistical dispersion intended to summarize the information contained in the [START_REF] Lorenz | Methods of measuring the concentration of wealth[END_REF] curve that represents the income distribution of a nation's residents; G may be viewed, though, more generally as a metric summarizing the dispersion of any smooth curve that starts at the origin and ends at the point [START_REF] Kondrashov | Log-normal Kalman filter for assimilating phase-space density data in the radiation belts[END_REF][START_REF] Kondrashov | Log-normal Kalman filter for assimilating phase-space density data in the radiation belts[END_REF] with respect to the diagonal of the corresponding square. In particular, we use G here to summarize into a single scalar the ROC curve, which ranges from 0 for random discrimination to 1 for perfect discrimination.

Figure 4a shows ROC curves obtained over the entire sample of n = 50 000 sequences: they correspond to G = 0.35 for the conventional method and to G = 0.82 for the DADA method, i.e. the overall performance gap is more than twofold. As expected, the performance of both methods is nil for λ 1 = 0 and it is very sensitive to the intensity of the forcing, cf. Fig. 4b.

Furthermore, the skill of the DADA method is boosted when decreasing the level of model error, cf. Fig. 4c; this is an expected result, since DA becomes more reliable when the model is more accurate, and when it is known to be so. Ultimately, under perfect model conditions, i.e. as σ Q → 0, DADA reaches perfect discriminative power, with G→ 1, no matter how small, but still positive, the forcing is; see Fig. 4d. On the other hand, the level of observational error σ R appears to have but a limited effect on DADA performance for the range of values considered, cf. Fig. 4e.

Finally, Fig. 4f shows that both methods perform better when the contrast between p 0 and p 1 is strong, but the latter does not influence the gap between the two methods, which remains nearly constant. This constant gap thus appears to quantify the additional power resulting from the extra discriminative features that the PDF f (y) is able to capture on top of those associated with the probability P (φ(y) ≥ u).

4 Discussion and conclusions [START_REF] Hannart | Counterfactual causality theory for the attribution of weather and climate-related events[END_REF] have relied on the causality theory of [START_REF] Pearl | Causality: Models, Reasoning and Inference[END_REF] to show that the ratio between the factual evidence f 1 (y) and the counterfactual evidence f 0 (y) is important in studying causal attribution of weather-and climate-related events. In this paper, we first described data assimilation (DA) methods and then demonstrated that they are well suited for deriving f 0 (y) and f 1 (y) from trajectories in the factual and the counterfactual worlds, respectively. Besides, these methods offer the key practical advantage of being already up-and-running in near real time at meteorological centers.

Combining these two sets of considerations, theoretical and practical, opens a novel route towards near real time, systematic causal attribution of weather-and climate-related events, thereby addressing a key challenge in the field of detection and attribution (D&A) at present [START_REF] Stott | Attribution of weather and climate-related events[END_REF].

Theoretical considerations

Implementing the DA for D&A (DADA) approach in the context of the L63 model in Section 3 allowed for a detailed step-by-step illustration of our methodological proposal. It also provided a basic test for an initial performance assessment, which showed an improved level of discriminating power with respect to the conventional approach outlined in Section 1. These results are promising, and their promise is easy to understand, given the fact that the DADA approach leverages the available information on the entire trajectory y, as opposed to the single specific feature 1 φ(y)≥u in the conventional approach.

It is important, though, to stress that the term "performance" here should be considered with caution: improving discriminatory performance may or may not be a desirable outcome, depending on the causal question being asked. [START_REF] Hannart | Counterfactual causality theory for the attribution of weather and climate-related events[END_REF] have shown that the causal question being formulated reflects the subjective interests of a particular class of end-users, and that the formulation itself may dramatically affect the answer.

For example, the question "did anthropogenic CO 2 emissions cause the heatwave observed over Argentina during January 2014? " has been traditionally treated by defining a "heatwave" in terms of a predefined temperature index reaching a predefined threshold, i.e., by a singular index exceeding a singular threshold. This class of questions matters for instance in the context of insurance disbursements, where a financial compensation may typically be triggered based on such an index exceedance. In this situation, the additional discriminatory power of DADA is meaningless because the DADA computation does not address the question at stake: there is simply no alternative to computing the probabilities p 0 and p 1 of the index exceeding the threshold.

However, if the question is formulated instead as "did anthropogenic CO 2 emissions cause the atmospheric conditions observed over Argentina during January 2014? " -i.e., without specifying which feature of the observed sequence is most important -then improving discrimination makes perfect sense and DADA becomes fully relevant. Furthermore, DADA is still fully relevant even if the question is formulated more specifically as "did anthropogenic CO 2 emissions cause the damages generated in Argentina by the atmospheric conditions of January 2014?," provided that is, that a model relating atmospheric observations to damages at every time step t along the trajectory of the physical model used in the assimilation is available and can be integrated into the observation operator H.

On the other hand, the results of Section 3 should also be considered with caution simply because the L63 testbed obviously differs in many respects from the real situation envisioned for future applications, both in terms of model dimension n and observation dimension d: in practice n will be very large and d n, while here we took d = n = 3.

In particular, choosing a highly idealized, climatological a priori distribution on the initial condition π(x 0 ) does not raise any difficulty under the tested conditions nor does it influence significantly the outcome of the procedure (not shown). The choice of π(x 0 ), however, may be an important problem in practice, when d n, and lead to potentially spurious results.

As a consequence, it may be both necessary and useful to further constrain the so-called background PDF π(x 0 ) by using the forecasts originating from τ previous assimilation cycles, thus following the ideas of laggedaveraged forecasting [START_REF] Hoffman | Lagged average forecasting, an alternative to Monte Carlo forecasting[END_REF][START_REF] Dalcher | The ensemble Kalman filter: theoretical formulation and practical implementation[END_REF]. The evidence thus obtained, though, will then also depend on previous observations over the "initialization" window [-τ, ..., -1] -i.e., it will no longer represent exclusively the desired evidence f (y). Besides, choosing τ optimally to constrain the initial background PDF in a satisfactory manner, while at the same time limiting the latter unwanted dependence on previous observations, is a challenging question that needs to be adressed.

More generally, the problem of evaluating the evidence f (y) is not new in the HMM and DA literature; see, for instance, [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF]; [START_REF] Hürzeler | Approximation and maximising the likelihood for a general state-space model[END_REF]; [START_REF] Pitt | Smooth particle filters for likelihood evaluation and maximisation[END_REF] and [START_REF] Kantas | An overview of sequential Monte Carlo methods for parameter estimation[END_REF]. Various algorithms are thus available to carry out this evaluation, depending on a number of key assumptions -such as lack of Gaussianity or linearity -and on the inferential setting chosen, e.g. particle filtering. These algorithms may provide accurate and effective solutions to the above problem, as well as improved alternatives to the Gaussian and linear approximation of Eq. ( 8), since the latter may not be sufficiently accurate for succesfully implementing the DADA approach under realistic conditions.

Practical considerations

While we have shown here that the proposal of using DADA for event attributions has intellectual merit, its main strength lies, in our view, in down-to-earth cost considerations. By design, the DADA approach allows one to piggyback at a low marginal cost on the large and powerful infrastructures already in place at several meteorological centers, in terms of both hardware and personnel. These centers are capable of processing massive amounts of observational data with highthroughput pipelines on the world's largest computational platforms, as opposed to requiring the design, set-up and maintenance of a new and large, D&A-specific infrastructure to collect observations and generateunder near real time constraints -the many model simulations required by the conventional approach recalled in Section 1.

Taking a step back, it is useful to examine our proposal within the wider context of the emergence of socalled climate services. It is widely recognized that extending the scope of activity of meteorological centers from being "monoline" weather forecasting providers to becoming "multiline" climate services providers -encompassing, for instance, weather forecasting and weather event attribution as two service lines among several others -is a relevant strategic option [START_REF] Hewitt | The Global Framework for Climate Services[END_REF]. Such a strategy may foster the timely and cost-efficient emergence of the latter services by building upon technological and infrastructure synergies with the former. For these reasons, our proposal is particularly relevant for, and could contribute to, the implementation of the strategic option just outlined.

This being said, DADA can very well serve as a method for near real time event attribution even for hypothetical climate services providers that focus uniquely or mainly on longer time scales, beyond a month, a season or a year. In such a context, DADA may allow for the assimilation of a broader range of observations, and in particular of ocean observations; it may, in fact, be important to include the latter in causal analysis when the event occurrence under scrutiny is defined over a sufficiently large time window. superscript f in the DA literature [START_REF] Ide | Unified notation for data assimilation: Operational, sequential and variational[END_REF]. It follows that this decomposition is tractable using a DA scheme that is able to estimate the conditional and forecast PDFs.

Next, let us apply the general Bayesian inference (10) to the case in which all the PDFs involved are Gaussian; this requires, in turn, that both the dynamics and observation models M and H be linear, and that the input statistics all be Gaussian. In this case, the Kalman filter allows for the exact computation of the PDFs mentioned in Eq. ( 10), which turn out to be Gaussian.

In the following, N (x, P) designates the Gaussian PDF of mean x and covariance matrix P. In this context, the analysis PDF at instant t is N (x a t , P a t ), where x a t and P a t are the analysis state and error covariance matrix at instant t. As a result of the linearity assumptions, the forecast PDF at instant t + 1 is given by a Gaussian distribution N (x f t+1 , P f t+1 ), where x f t+1 and P f t+1 are the forecast state and error covariance matrix at instant t + 1. Further, the integration on x t+1 in Eq. ( 10) can readily be performed under these circumstances, with the outcome that p(y t+1 |y t ) is distributed as N (H x f t+1 , R + H P f t+1 H ). The desired model evidence f (y) can then be computed by recursion on successive time steps as:

f (y) = p(y 0 ) T t=1 (2π) -d 2 |Σ t | -1 2 × exp - 1 2 (y t -Hx f t ) Σ -1 t (y t -Hx f t ) ; (11)
here p(y 0 ) represents the prior PDF of the initial state, Σ t = R + HP f t H , and This expression coincides with Eq. ( 8) and can be evaluated with the help of any DA method that yields the forecast states and forecast error covariance matrices, such as the KF or the EnKF. Note that the traditional standard Kalman smoother would give the same result as the KF, since they share the same forecasts.

Finally, Eqs. ( 10) and ( 11) above show that the likelihood f (y) may be obtained as a by-product of the inference on the state vector x, which usually is the main purpose in numerical weather prediction. This idea may actually be highlighted in even greater generality by considering the equality:

f (y) = p(y | x)p(x) p(x | y) . (12)
While Eq. ( 12) is a direct consequence of Bayes theorem, it also illustrates a point that is arguably not so intuitive. The likelihood f (y) is obtained here as the ratio of two quantities: a numerator p(y | x)p(x) that is a model premise inherently postulated by Eqs. ( 5) and ( 4), and a denominator p(x | y) that may be viewed as the end result of the primary inference on x. In other words, estimating f (y) requires only a straightforward division, provided x has been previously inferred. Equation ( 12) thus expresses with great clarity and simplicity a fundamental idea buttressing our proposal, as it provides a general theoretical justification for the suggestion of deriving the likelihood from an inferential treatment that focuses on x. To put it succintly, this equation basically says, "He who can do more can do less." In the context of DA, whose end purpose is to infer the state vector x out of an observation y -i.e., the more part -it is possible to obtain the likelihood as a by-product thereof -i.e., the less part -and thus almost for free.

Fig. 1 Fig. 2

 12 Fig. 1 Illustration of the conventional D&A approach as applied to a univariate AR(1) process. (a) Observed time series (first component Y 1 , dotted line) and daily average φ(Y) (heavy solid line). (b) Threshold level (vertical axis) as a function of the return period (horizontal axis): simulated values (crosses); fit based on the Generalized Pareto distribution (GPD, heavy dark-blue line); uncertainty range at the 95% level (light blue area); and threshold value u = 3.1 (light solid black line). (c) Estimated value of P = P (φ(Y) ≥ u) (heavy dark-blue line) using a GPD fit as a function of the sample size n (horizontal axis); uncertainty range (light blue area); and true value P = 0.01 (light solid black line). (d) Computational time on a desktop computer (seconds, vertical axis) as a function of sample size n (horizontal axis) required by the conventional method (dark blue line) and the DADA method (solid red line); the latter method is explained in Sections 2b and 3 below.

Fig. 3

 3 Fig. 3 Sample trajectories from data assimilation (DA) in our modified L63 model. (a) True trajectory (black solid line) and the two trajectories reconstructed by DA in the factual (i = 1) and counterfactual (i = 0) worlds (red and green solid lines), respectively, over a long sequence, T = 400; the values of λ 1 and θ 1 here are the same as in Fig. 2, and the assimilated observations are shown as gray dots. (b) Same as panel (a) but zoomed over a short sequence, T = 20. (c) Logarithm of the cumulative evidences f 1 (y) and f 0 (y) (red and green lines, respectively) computed over the window [0, t ≤ T ]; gray bars indicate the instantaneous differences between f 1 (y t ) and f 0 (y t ). (d) PN computed over the window [0, t].

Fig. 4

 4 Fig. 4 Performance of the DADA and conventional methods (red vs. blue solid lines, respectively). (a) Receiver operating characteristic (ROC) curve: true positive rate as a function of false positive rate, when varying the cut-off level u, as obtained from the entire sample of n = 50 000 sequences; see text for details.. (b) Gini index G as a function of forcing intensity λ 1 . (c) Same as (b) for several values of σ Q and for DADA only, with the black arrow indicating the direction of growing σ Q . (d) Same as (b) but as a function of model error amplitude σ Q . (e) Same as (b) but as a function of observational error amplitude σ R . (f) Same as (b) as a function of the logarithmic contrast between the conventional probabilities log p 1 /p 0 .

Acknowledgements This work has been supported by grant DADA from the Agence Nationale de la Recherche (ANR, France: AH and all co-authors) and by the Multi-University Research Initiative (MURI) grant N00014-12-1-0911 from the the U.S. Office of Naval Research (MG).

Appendix: Derivation of the model evidence In this appendix, we outline the derivation of model evidence within a general Bayesian framework, and we apply the latter to the narrower KF context to obtain Eq. ( 8). Consider two consecutive cycles of a DA run, the first with state vector x t and observation vector y t at instant t and the subsequent one with state vector x t+1 and observation vector y t+1 at instant t + 1. We plan to find a tractable expression for the model evidence p(y t , y t+1 ).

The model evidence provided by the full sequence of observations y = (y 0 , ..., y T ) will be inferred by recursion, using the results of this two-observation setting. In order to decouple the two cycles, one first has to spell out the Bayesian inference p(y t , y t+1 ) = p(y t )p(y t+1 |y t ). We look for a tractable expression for p(y t+1 |y t ) by further introducing the states x t+1 and x t as intermediate random variables:

where p(y t+1 |x t+1 ) is the likelihood of the observation vector y t+1 conditional on the state vector x t+1 and it is known from Eq. ( 5).

The conditional PDF p(x t |y t ) of x t on y t at instant t -which appears on the right-hand side of the above equation -is referred to as the analysis PDF in the DA literature, where it is denoted by a superscript a [START_REF] Ide | Unified notation for data assimilation: Operational, sequential and variational[END_REF], and it constitutes the main DA output. The integral xt p(x t+1 |x t )p(x t |y t ) dx t = p(x t+1 |y t ), in which p(x t+1 |x t ) is known from the model dynamics given by Eq. ( 4), propagates this analysis PDF further in time, to instant t + 1. Hence, the result of this integration coincides with the forecast PDF, denoted by