A D El Maliani

M El Hassouni

Y Berthoumieu

D Aboutajdine

Color texture classification method based on a statistical multi-model and geodesic distance

Keywords: Color texture, wavelet representation, Rao geodesic distance, Bayesian classification

HAL is a

Introduction

For color texture analysis, the question of the color space selection to characterize the visual content may be raised. RGB space is the most natural since it represents the acquisition-reference space. However, R, G and B channels are known to exhibit simultaneously high statistical dependencies relatively to: 1) the intra-channel or spatial one, 2) the inter-band for scale-space representation one and 3) the inter-channel or multicomponent one. In recent works, some authors have proposed various stochastic models to characterize the inter-channel dependencies alone. From scale-space representation such as 2-D wavelet or others, previous approaches consist in fitting the histograms of the wavelet subbands with a given parametric probability density function (pdf) of the trivariate color vector samples for each level of the decomposition. Recent works are based on multivariate generalized Gaussian distribution in [Verdoolaege et al. 2012] or copulas based multivariate Weibull distribution in [Kwitt et al. 2011]. The parameters of the pdf are the covariance matrix and different auxiliary parameters such as the shape parameter characterizing the non-Gaussiannity. Extending these models to take into account jointly intra-band, inter-band and inter-channel dependencies, leads to a high-dimensional model in terms of number of parameters. This is a main drawback of these proposals since, in practice, for a fixed data size, dealing with high-dimensional model may be a critical issue in terms of parameter estimation performance. Considering color space, another alternative to the RGB space is to consider the family of luminance-chrominance (LC) color spaces [DeYoe et al. 1996, Mojsilovic et al. 2000]. This family includes perceptual (HSV, HSI and HSB) and perceptually uniform (Lab, Yuv, YCbCr) color spaces [Sarifuddin et al. 2005]. An extended comparison of those spaces is given by [Qazi et al. 2011]. These authors note that the luminance channel is uncorrelated with the chrominance channels. From this independence in LC spaces, they proposed a model taking into account separately intra and inter dependencies, i.e. multi-model, with few parameters compared to the required number by using the RGB space. Qazi et al. proposed to consider a multi-model based on a 2-D complex quarter plan autoregressive (2-D QP AR) model for luminance as well as for chrominance information. This may lead to a lack of flexibility because the luminance and chrominance are decorrelated and thus, have distributions of different natures.

In this letter, considering the LC color spaces, we propose a multi-modeling approach based on multivariate stochastic modeling. The contribution is twofold:

• First, we propose a multi-model based on two non-Gaussian sub-models, i.e. M L and M Cr , devoted to describe respectively 1) the intra-channel dependency by considering the spatial structure and inter-band dependence of the luminance channel, and 2) the inter-channel dependency by considering a bivariate pdf for the two chrominance channels. Both M L and M Cr are based on copula as a powerful tool to construct non-Gaussian joint models that fit data with even different natures of margins. Linear and/or circular models are derived depending on the component natures (ex: the Hue in HSV color space is a circular channel).

• Second, for the similarity measurement used to compare two texture samples, we derive a closed form expression for the Rao geodesic distance for copula based multivariate models. Contrary to the conventional Kullbackleibler (KL) divergence, the Rao geodesic metric has the advantage to satisfy distance properties such as symmetry and triangular inequality [Atkinson et al. 1981].

The use of different models for characterizing the wavelet coefficients has been recently investigated by [Rakvongthai et al. 2013], for texture retrieval in a noisy environment. Authors proposed different univariate statistical models for the real/imaginary, magnitude and phase parts of complex subband coefficients.

The parameter vector of the studied texture is then formed by the concatenation of parameters estimated from each of these models. Also, the similarity measure has been expressed by a formula combining both magnitude and phase KL divergences using a mixing parameter. While the cited work did not consider the color information at all, it should be emphasized that the present work considers the color and texture information jointly for image characterization.

Concerning the similarity measure, the originality of our method is the use of Rao geodesic distance which represents a good alternative to the KL divergence since it satisfies the distance properties. For this, a closed-form expression of Rao geodesic distance between copula based distributions is developed.

Note also that this work stands out from our previous work [Maliani et al. 2012] which presents a preliminary study of the multi-modeling concept that considered only the special case of color representation, the HSV color space. Moreover, the paper did not discuss the dimensionality reduction issue when considering the three dependencies (intra-band, inter-band and inter-channel). Also, in the classification step, we proposed the K-Nearest Neighborhood classifier taking the weighted sum of the two similarity measures between the luminance and chrominance sub-models. The weights have been tuned experimentally which is the same as for [Rakvongthai et al. 2013], which led us to test exhaustively on several scales of weights to get the best classification rate. The multi-modeling proposed in this letter, helps in dimensionality reduction of the extracted features to characterize the textural content. According to several dependence experiments presented in section 2, the spatial structure of the texture and the inter-band dependence is only considered for the luminance component and the inter color-channel dependency consists of a bivariate model. The performance of the multi-model is tested using a Bayesian classifier based on the softmax principal [Bishop 2006] thanks to its suitability to combine several similarity measures through a posteriori probability fusion. The experimental protocol adopted here, aims to find the optimal configuration of our multi-model. This starts with the choice of statistical distribution related to color components in the field of dual-tree complex wavelet CT-DWT [Kingsbury 1998]. We proposed to build joint models by testing different combinations of Gamma , Weibull and VonMises distributions for both luminance and chrominance components. First comparisons were considered to show the relevance of each model for the color spaces. Then, we compared our best selected models with the state of the art methods in terms of classification rate, the dimensionality reduction and the computation time.

This paper is organized as follows. Section 2 deals with the proposed multimodeling approach for textures in LC color spaces. In Section 3, we present the classification framework, and the proposed similarity measurement. Then, we show experimental results in Section 4. Finally, we give concluding remarks in Section 5.

Color texture modeling in LC color spaces

Dependence characterization

Let us consider a stationary color texture x T represented in a given LC color space, and let l be the luminance component, cr 1 /cr 2 respectively the first and the second chrominance component of x T . Each color component is first decomposed into a set of wavelet subbands. In wavelet domain texture characterization, it is common to model the statistical distribution of wavelet coefficients in different subbands, which corresponds to the detail subbands. The lowpass wavelet coefficients are excluded ( [Mallat 1999], [START_REF] Boubchir | [END_REF]). The subbands are respectively l (k) , cr

(k) 1 , cr (k) 2 , k = 1, • • • , K.
The statistical distribution of wavelet detail coefficients of natural textures has been extensively studied ( [Simoncelli 1997], [START_REF] Boubchir | [END_REF]], [START_REF] Lasmar | [END_REF]], [Verdoolaege et al. 2012], [Kwitt et al. 2011]). It has been shown that the empirical distribution is generally heavy-tailed with a sharp peak at zero. This leptokurtic behavior leads to consider non-Gaussian marginal or joint density.

Contrary to the marginal approach, the joint modeling needs to define the dependence which exists between subbands coefficients. Intra-band, inter-band and inter-channel dependencies have to be considered. Let us detail each of these dependencies for the three color components.

Intra-band dependency

Intra-band dependence refers to the spatial neighborhood dependence between coefficients of the same wavelet subband. It corresponds to the spatial structure of the studied subband. Let s be the location of a reference coefficient from the k th subband of a given color channel. If we consider the luminance channel for example, let m be the size of the chosen neighborhood. The collection of the reference coefficient neighbors is specified by l

(k) s = l (k) 1,s , • • • , l (k) m,s .
Then, considering all the subbands coefficient neighbors, the dataset to be observed is:

l intra = l (k) 1 , • • • , l (k) s , • • • , l (k) P T (1)
l intra is a P by m matrix, where P is the size of the wavelet subband. If we consider, the first chrominance channels, intra-band dependency for a given subband k is studied by observing the following dataset:

cr 1intra = cr 1 (k) 1 , • • • , cr 1 (k) s , • • • , cr 1 (k) P T (2)
For the second chrominance channel, we have:

cr 2intra = cr 2 (k) 1 , • • • , cr 2 (k) s , • • • , cr 2 (k) P T (3)

Inter-band dependency

To model the inter-band dependence, we consider subbands in different orientations for the same color channel. Thus for luminance inter-band we observe a dataset represented by a P by K matrix, where K is the total number of subbands:

l interbd = l (1) • • • l (K) (4)
for chrominance inter-band:

cr 1 interbd = cr (1) 1 • • • cr (K) 1
(5)

cr 2 interbd = cr (1) 2 • • • cr (K) 2 (6)
we note that l (k) , cr

D interch = l (k) cr (k) 1 cr (k) 2 (7)

Full dependence

Considering all the aforementioned dependencies leads to a full dataset that considers intra-band, inter-band and inter-channel dependencies jointly. It is represented by the following matrix of size P by (3 × K × m):

D = l (1) 1 , • • • , l (1)
P , • • • , l (K) 1 , • • • , l (K) P , cr 1 (1) 1 , • • • , cr 1 (1) P , • • • , cr 1 (K) 1 , • • • , cr 1 (K) P , cr 2 (1) 1 , • • • , cr 2 (1) P , • • • , cr 2 (K) 1 , (8)
• • • , cr 2 (K) P T
Consideration of this matrix means, first, that all the three components carry the spatial structure information (intra-band). Second, that dependence exists between subbands of all orientations for each channel (inter-band), and third, that the representation exhibits dependence among the three color channels (inter-channel).

Multi-model approach

Experimental analysis of dependence

Contrary to the RGB representation, LC color spaces offer separability between luminance and chrominance. At first glance, this indicates that the interchannel dependence have not to be considered between luminance and chrominance. A more precise study of dependence is done according to the correlation coefficients. Table 1 shows degrees of dependence when textures are represented in Lab color space. This is obtained by averaging the Pearson, Kendall, Spearman and the mutual information coefficients on all subbands of 40 color textures from the Vistex database [Vistex]. We recall that mutual information is a measure of the statistical dependence between two variables. If we suppose that X and Y are two random variables, f (x) and f (y) the marginal probability distribution functions of X and Y respectively, and f (x, y) the joint probability distribution function of X and Y . The mutual information between X and Y is then:

I(X, Y ) = x,y f (x, y)log f (x, y) f (x)f (y) (9)
We deduce from these results that intra-band dependence have to be consid- ered only for luminance, and that inter-channel dependence exists only between chrominance components. They also show that we can restrict the study of inter-band dependence to the luminance component.

Multi-model formulation

Our contribution here is to exploit the independence between luminance and chrominance. So we propose a multi-model based characterization that considers the different natures of these information ( [Qazi et al. 2011]). Two joint and independent models M L and M Cr are used to respectively characterize the luminance information with spatial and inter-band dependence considering observations of a vector l s (see subsection 2.4.1) and the chrominance one devoted to inter-channel dependence by observing realizations of a vector cr s (see subsection 2.4.2):

p(x T , θ) = p(x T |{M L , M Cr }) ∝ p(l s , θ L ) K k=1 p(cr (k) s |θ cr ) ( 10
)
where θ L and θ cr are the vectors of parameters of luminance and chrominance models, respectively. The multi-modeling permits to reduce the dimension of the observed dependence dataset D (subsection 2.1).

Copula based multivariate models

Different multivariate models can be used for defining M L and M Cr , such as the multivariate generalized Gaussian (MGGD) [Verdoolaege et al. 2012], the multivariate Bessel K forms [START_REF] Boubchir | [END_REF], the spherically invariant random vector (SIRV) [START_REF] Lasmar | [END_REF]]. Those models have been proposed to characterize different kinds of dependencies such as the intra-band in the case of SIRV and the inter-channel in the case of MGGD. Each of them is dedicated to a specific dependency rule, either for intra-band or inter-channel but never taking account of both together. Obviously, extending each model simultaneously to the both dependencies is possible but the equivalent model exhibits the main drawback the common definition of the leptokurtic and platykurtic for the whole set components of the multivariate model. If we suppose to have a

set of dependent variables T = [t 1 , • • • , t n ],
with the MGGD or SIRV multivariate modeling, each marginal t i has finally a similar shape value. This behavior is confirmed in the study of [Rangaswamy et al. 2002], the characterization of one marginal is sufficient to determine the shape of the whole distribution of a SIRV. Table 2 show averages of scale and shape parameters supposing that subbands of the three components are described by the same marginal, Weibull or Gamma for example. It is clear from the table that the marginal modeling of luminance and chrominance exhibits different shape and scale parameters in the case of Lab color space, contrary to RGB. In the HSV color space, we observe that parameters are different even between the chrominance components H and S. The copula theory allows the joint modeling by separating the dependence structure from the marginal modeling. It is thus easy to merge different marginal laws into the same multivariate pdf. [Nelsen 2006] defined the copula as a multivariate uniform distribution used to construct multivariate models. Given that F is a d-dimensional cumulative distribution function (cdf) with continuous marginals F 1 , ..., F d , the Sklar's theorem [Sklar et al. 1959] shows that there exist a unique copula C such that:

F (x 1 , • • • , x d ) = C (F 1 (x 1 ), • • • , F d (x d )) (11)
Further, if C is continuous and differentiable, the copula density (c) is given by:

c(u 1 , ..., u d ) = ∂ d C(u 1 , ..., u d ) ∂u 1 ...∂u d (12)
The joint pdf is uniquely deduced from the marginals and the copula density as follows:

f (x 1 , ..., x d ) = c (F 1 (x 1 ), ..., F d (x d )) d i=1 f i (x i ) ( 13
)
where

f i , i = 1, • • • , d, represent the marginal pdfs that can follow different laws.
Estimation of the copula parameters is done using the canonical maximum likelihood (CML) [START_REF] Durrleman | [END_REF]] method. This is a pseudo parametric estimation method contrary to the fully parametric method, termed the inference from marginals (IFM) [Joe. 1997]. The CML method consists in transforming the data (x 1 , ..., x d ) into uniform data (û 1 , ..., ûd ) using the empirical distributions (empirical cdfs) and then estimates the copula parameters considering all the observations as:

Σ = argmax Σ n i=1 log c (û i1 , ..., ûid ; Σ) (14)
where n denotes the number of observations, and Σ represents the correlation matrix.

For marginal parameters, the estimation is done using the conventional maximum likelihood method. Experiments reveal that the CML estimator is the best since the method does not rely on parametric assumption about marginals.

This means that less pertinent choice of marginals do not affect estimation of the copula parameters.

For texture modeling, we repose on copulas to construct the models M L and M Cr . We select the Gaussian copula for which the density is defined by:

c φ (u 1 , • • • , u d ) = 1 |Σ| 1/2 exp - 1 2 ϑ T (Σ -1 -I)ϑ (15) with ϑ T = φ -1 (u 1 ), • • • , φ -1 (u d ) ,
φ denotes the cdf of the normalized Gaussian density, Σ the correlation matrix, and I the d-dimensional identity matrix.

The choice of the Gaussian copula is justified by:

• the dependence representation, which is based on the well known correlation coefficient.

• the existence of closed form of Rao distance for the Gaussian copula which allow us to derive a closed form of Rao distance for the joint model as a similarity measure (subsection 3.1).

M L and M Cr models

Considering the full dependence dataset D (see sub-section 2.1.4) leads to high parametrization and thus to very cumbersome estimation of parameters.

To alleviate this complexity, we use the proposition made by [Qazi et al. 2011] about dependencies between coefficient subbands when color textures are represented in LC color spaces. The proposition consists in a multi-modeling approach to reduce the dimensionality of the observed dataset. Considering all positions at the luminance subbands, the dataset to be modeled is a P by d matrix of realizations of the vector l:

L =               l 1 l 2 . . . l s . . . l P              
where P is the size of a wavelet subband. The observation matrix is:

L =            l 1,1 • • • l m,1 • • • l d,1 . . . . . . . . . l 1,s • • • l m,s • • • l d,s . . . . . . . . . l 1,P • • • l m,P • • • l d,P           
Based on the Gaussian copula, the pdf of the model M L is defined by:

f M L (l s ; θ) = 1 |Σ L | 1/2 exp - 1 2 ϑ T s (Σ -1 L -I)ϑ s d i=1 f s (l i,s ) (16)
where θ = (w 1 , w 2 , ..., w d ; Σ L ) is the vector of hyperparameters,

w s = {w sj } j∈ℵ , s = 1, • • • , d
is the vector of parameters for the marginal f s and Σ L is the covariance matrix of the Gaussian copula.

In the case of Gaussian copula, the CML estimator of covariance matrix is the sample correlation matrix of Gaussian observations ϑ 1 , • • • , ϑ P :

ΣL = 1 P P s=1 ϑ s ϑ T s ( 17
)
with

ϑ s = (ϑ 1,s , • • • , ϑ d,s
) is obtained by transforming the observation

l s = (l 1,s , • • • , l d,s ) from the observation matrix L by ϑ i,s = φ -1 (F i (l i,s )), i = 1, • • • , d
and F s is the empirical marginal cdf.

2.4.2. Circular/Linear bivariate model for chrominance (M Cr )

We consider the inter-channel dependence only between chrominance component subband represented by the vector cr = cr

(k) 1 , cr (k) 2
. Thus, considering all coefficients of the chrominance subbands, the dataset to be modeled by M Cr is a P by 2 matrix:

cr =         cr 1,1 cr 2,1 . . . . . . . . . . . . cr 1,P cr 1,P        
M Cr is a bivariate model defined reposing on the Gaussian copula:

f M Cr (cr (k) s ; θ) = 1 |Σ Cr | 1/2 exp - 1 2 ϑ T s (Σ -1 Cr -I)ϑ s 2 i=1 f i (cr i ) ( 18
)
where θ = (w 1 , w 2 ; Σ Cr ), w i = {w ij } j∈ℵ . The index i ∈ {1, 2} is the vector of parameters for the marginal f i , and Σ Cr the covariance matrix. The CML estimator of the covariance matrix is as follows:

ΣCr = 1 P P s=1 ϑ s ϑ T s ( 19
)
with ϑ s = (ϑ 1,s , ϑ 2,s ) is obtained by transforming the observation cr s = (cr 1,s , cr 2,s )

from the observation matrix cr by ϑ i,s = φ -1 (F i (cr i,s )), i = 1, 2 and F s is the empirical marginal cdf.

The model is called circular/linear because it offers the possibility of using different marginals f i to come up with a bivariate circular/linear, circular/circular or linear/linear model for chrominance channels depending on the nature of the used color space.

Texture joint modeling

Taking into account the previous sub-models, respectively M L and M Cr termed probabilistic because they describe probability distributions over subband coefficients of the texture, let us consider a stationary observed texture, we have:

p(x T |{M L , M Cr }) ∝ p(l s |w s , Σ L ) K k=1 p(cr (k) s |w (k) cr1 , w (k) cr2 , Σ (k) cr ) ( 20
)
which can be written in a more compact parametric form as:

p (x T |{M L , M Cr }) ∝ p(x T |θ s ) (21) with θ s = w s , Σ L , {w (k) cr1 , w (k) cr2 Σ (k) cr } k=1,•••,K

Classification

In the framework of stochastic parametric classifier, we consider supervised representing the class and the other the texture query x T , is expected. Thus, the class assignment is obtained from a similarity measure. To do this, the Softmax approach or normalized exponential function [Bishop 2006] is selected. in the seminal theory of the information geometry. This consists in considering that the parametric model of finite dimension forms a sub-manifold characterized by a smooth Riemannian space. Some of these measures are divergences such as Kullback-Leibler. A divergence is a non-symmetric measure and can be not respecting the triangle inequality. Some others such as Bhattacharya, Hellinger or others are really some distances. Considering the properties of the Riemmanian manifold, the well-founded distance associated to this type of manifold is the geodesic or Rao distance [Atkinson et al. 1981]. The main drawback of the geodesic distance is that a closed form for probabilistic model is rarely available. This is the reason for which the Kullback-Leibler divergence is usualy used in place to the geodesic distance. However, the true distance according to the Riemannian space is the geodesic distance.

In the next subsection, we derive for our context a closed form for the geodesic distance for the Copula modeling.

Proposed Rao distance as similarity measure

For a probability density f (x; θ) with θ = (θ 1 , • • • , θ N ) a vector of parameters, the Rao distance is a Riemannian metric defined by the fisher information matrix as:

ds 2 = N p,q=1
g pq (θ)dθ i dθ j (24)

where g pq represent the Fisher matrix elements:

g pq (θ) = -E ∂ 2 ∂θ p θ q log f (x; θ) (25)
Let f be a Gaussian copula based joint pdf (M L or M Cr ), θ = {w i , Σ} is the vector of hyper-parameters of f , w i = {w ij } j∈ℵ is the vector of parameters for the marginal f i , i = 1, • • • , d (for example if the marginal is a Weibull we have w i = {α i , β i } representing shape and scale parameters), d is the length of the vector of one observation and Σ represents the covariance matrix of the Gaussian copula. By applying equations 13 and 15 in equation 25 we have:

g pq (θ) = -E ∂ 2 ∂θ p θ q log c Φ (u; Σ) + log d i=1 f i (x i ; w i ) ( 26
)
where u = (u 1 , ..., u d ) represents the vector of uniform coefficients that are obtained by transforming the wavelet coefficients (x 1 , ..., x d ) using their cdfs (u i = F -1 (x i , w i )) or using their empirical cdfs (see MATLAB's routine "ecdf"

) (u i = eF -1 (x i ))
. In this latter case, no assumptions are made on the parametric form of the marginal distributions. The copula parameters are then independent from any marginal parameter w ij :

∂ ∂w ij c Φ (u; Σ) = 0 (27) so then, g wij Σ (θ) = g Σwij (θ) = 0, i = 1 • • • d, j ∈ ℵ Thus, g ΣΣ (θ) = -E ∂ 2 ∂Σ∂Σ log c Φ (u; Σ) ( 28
)
and

g µν (θ) = -E ∂ 2 ∂µ∂ν log d i=1 f i (x i ; w i ) (29) with µ, ν ∈ {w ij }, i = 1 • • • d, j ∈ ℵ g µν (θ) = -E ∂ 2 ∂µ∂ν d i=1 log f i (x i ; w i ) (30) = -E d i=1 ∂ 2 ∂µ∂ν log f i (x i ; w i ) (31) = d i=1 -E ∂ 2 ∂µ∂ν log f i (x i ; w i ) (32)
Thus, from (12):

ds 2 = g ΣΣ dΣdΣ + d i=1 µ,ν g µν μ νdt (33) = ds 2 Gauss + d i=1 ds 2 Margins (34)
Hence, given two copula based probability distributions f (1) (x; θ 1 ) and f (2) (x; θ 2 ) with θ 1 = (w

(1) 1 , • • • , w (1)
d , Σ 1 ) and θ 2 = (w

(2) 1 , • • • , w (
L 2 = L 2 Gauss (f (1) (x; Σ 1 )||f (2) (x; Σ 2 ))+ d i=1 L 2 Margins (f (1) i (x; w (1) i )||f (2) i (x; w (2) i )) (36)
that is:

L 2 = 1 2 d i=1 (ln λ i ) 2 + d i=1 L 2 Margins (f (1) i (x; w (1) i )||f (2) i (x; w (2) i )) (37)
where λ i , i = 1, ..., d represents the eigenvalues of Σ -1 1 Σ 2 .

EXPERIMENTAL RESULTS

The experimental section is concerned with five major issues:

1. Choice of the multi-model for color texture: for this, we investigate a variety of models for different color spaces.

2. Choice of the LC color space: we compare performances between two subfamilies of LC color spaces, namely the perceptually uniform color spaces (represented by Lab), and the perceptual color spaces (represented by HSV).

3. Comparison with existing approaches: we evaluate performances of our approach in comparison with existing uni-modeling approaches in LC color spaces [Qazi et al. 2010, Qazi et al. 2011] and in RGB color space [Verdoolaege et al. 2012, Kwitt et al. 2011].

4. Runtime: we address the runtime issue by comparing the proposed Rao distance based similarity measure with the Monte-carlo based Kullbackleibler divergence between copula based pdfs.

5. Curse of dimensionality: we compare our approach with the full feature vector based approach (subsection 2.1).

Experiments versus data diversity

In order to conduct representative experiments, we use different databases and various configurations of the experimental protocol leading to a large view of the proposed multi-model performance. The conventional Vistex database [Verdoolaege et al. 2012, Kwitt et al. 2011, Do et al. 2002] is used considering two sizes of samples respectively 32 × 32 and 128 × 128. Experiments are also conducted on the Outex, ALOT and Stex databases which are more challenging color texture database, since the color and texture information are not easily distinguishable. We consider the following scenarios:

• First scenario, i.e. DB1, addressing a set of 24 textured images of size 512×512 from the Vistex database [Vistex]. The protocol follows the work of [START_REF] Permuter | [END_REF]] and [Qazi et al. 2010]. Each image was divided into subimages of size 32 × 32 pixels resulting of a large database of 6144 small textures. We consider 96 image from the resulting 256 subimages as the training set, while the remaining 160 subimages are considered as the test set. This dataset is used in order to evaluate robustness of our model even if very local spatial structures are considered.

• Second scenario, i.e. DB2, addressing a set of 54 textured images of size 512 × 512 from the Vistex database. The protocol follows the work of [START_REF] Do | [END_REF]. Each image was divided into subimages of size 128 × 128 pixels. This dataset is available on the Outex web site [Ojala et al., Outex] as test suit Contrib TC 00006. For each texture, subimages are considered to form a checkerboard. The white half of subimages is then considered

• True positives (TP): refers to the number of subimages that have been successfully returned by the classifier as belonging to a given class.

• False positives (FP): is the number of subimages that have been wrongly classified as belonging to a given class.

• False negatives (FN): is the number of subimages that are classified as not belonging to a given class while, in fact, they belong to this class of textures.

• True negatives (TN): is the number of subimages that are correctly classified as not belonging to a given class.

The percentage of classification (or Sensitivity) is the proportion of textures which were well labeled by the classifier. It represents the ability of the classifier in correct positive assignment of a texture to its real class. A 100% percentage means that the classifier makes no false negative labeling (FN = 0). The Percentage classification is given as follows:

Sn = TP TP + FN × 100% (38)
However, theoretically, a system that always returns positive assignment (even wrongly) will achieve 100% percentage classification. The problem with percentage classification is that it tells nothing about textures that were wrongly labeled as belonging to a given class (FP assignments). For this, we take into account also the precision to evaluate the classifier performance.

The precision criterion indicates the ability of the classifier to detect true class membership. The precision is given by:

Sp = TP TP + FP (39)
A classifier with very high precision makes no mistakes in assigning a texture to its class (FP = 0).

Results and discussions

One of the advantages of the proposed multi-modeling approach is the flexibility in term of modeling textures from different LC color spaces. The joint

f (x, θ) = 1 |Σ| 1/2 exp[-1 2 ϑ T (Σ -1 -I)ϑ] × ( β -α Γ(α) ) d d i=1 x α-1 i exp -d i=1 ( x i β ), θ = (α, β, Σ)

CopMWbl

Gaussian copula based multivariate Weibull (all marginals are Weibull)

f (x, θ) = 1 |Σ| 1/2 exp[-1 2 ϑ T (Σ -1 -I)ϑ] × ( τ λ ) d d i=1 x τ -1 i exp -d i=1 ( x i λ ) τ , θ = (τ, λ, Σ)
{Gamma,Weibull} Gaussian copula based bivariate model with Gamma as the fisrt marginal and Weibull as the second marginal

f (x, θ) = 1 |Σ| 1/2 exp[-1 2 ϑ T (Σ -1 - I)ϑ] × τ β -α x α-1 2 λΓ(α) ( x 1 λ ) τ -1 × exp -( x 1 λ ) τ -x 2 β , θ = (α, β, τ, λ, Σ)
{vonMises,Weibull} Gaussian copula based circular/linear bivariate model with vonMises as the fisrt marginal and Weibull as the second marginal

f (x, θ) = 1 |Σ| 1/2 exp[-1 2 ϑ T (Σ -1 -I)ϑ] × τ 2πλI 0 (ν) ( x 2 λ ) τ -1 exp[νcos(x1- µ) -( x 2 λ ) τ ], θ = (µ, ν, τ, λ, Σ)
{vonMises,Gamma}Gaussian copula based circular/linear bivariate model with vonMises as the fisrt marginal and Gamma as the second marginal when chrominance components are defined in a cylindrical space, the hue represents a circular component and the saturation is the axe of the cylinder. We can then benefit from the advantage of copulas in merging different marginals in the same joint model to construct a circular/linear bivariate model for chrominance such as {vonMises,Weibull} or {vonMises,Gamma} (see Table 3).

f (x, θ) = 1 |Σ| 1/2 exp[-1 2 ϑ T (Σ -1 -I)ϑ] × β -α 2πΓ(α)I 0 (ν) ( x 2 β ) α-1 exp[νcos(x1- µ) -( x 2 β ) α ], θ = (µ, ν, α, β, Σ)
Regarding the first part of our experiments, Table 4 presents color texture classification performance in Lab and HSV color spaces respectively, using different multi-models on DB1, DB2, DB3, DB4 and DB5. We note that:

• Lab1: refers to a multi-model that considers CopMGam (Table 3) for the luminance model M L and CopMWbl for the chrominance model M Cr in Lab color space. We remark that in Lab, and for all the five databases, best performance is achieved using the multi-model Lab1. This indicates that the spatial structure along with the inter-band dependency are better characterized using the Gaussian copula in conjunction with Gamma marginals (CopMGam). This also indicates that the dependency between "a" and "b" chrominance components is better characterized using the Gaussian copula in conjunction with Weibull marginals (bivariate CopMWbl). For HSV, the multi-model HSV1 leads to the best performance in all databases. This is due to the ability of CopMGam to model spatial structure and inter-band dependence for luminance, and the circular/linear bivariate {vonMises,Weibull} to characterize the chrominance dependency.

We next would like to compare between the perceptually uniform and perceptual subfamilies of LC color spaces, via their respective representatives Lab and HSV spaces. Considering the best multi-models in each case (Lab1 for Lab and HSV1 for HSV), we clearly observe from Table 4 that Lab representation outperforms the HSV one. This is a fair comparison even the models M L and M Cr are not the same in multi-models Lab1 and HSV1, since these latters are the best relatively the nature of the color space. Then, if Lab1 outperforms HSV1, we can conclude that in best modeling case the Lab representation is better than the HSV one. These results are explained by the better luminancechrominance representation of Lab. In other words, we can say that Lab offers more decorrelation property between luminance and chrominance, which is the central interest of our approach. We do not ignore that percentage classification and precision rates are lower for large databases DB3, DB4 and DB5. Regarding the effectiveness of our approach in comparison with existing ones, i.e uni-modeling in LC and RGB color spaces, Table 5 shows average percentage classification Sn of textures of DB1, DB2, DB3, DB4 and DB5 using our best multi-models (Lab1 and HSV1) against the four following approaches:

• [Qazi et al. 2010], in Lab: Linear prediction models based characterization in the Lab color space,

• [Qazi et al. 2011], in IHLS: Linear prediction models based characterization in the improved hue luminance saturation (IHLS ) color space,

• [Verdoolaege et al. 2012], in RGB: Dependence between R, G and B color components using the MGGD model,

• [Kwitt et al. 2011], in RGB: Dependence between R, G and B using the t-copula based multivariate Weibull model, We observe from Table 5 that our Lab1 proposed model leads to better performance when compared with the linear prediction modeling [Qazi et al. 2010] and [Qazi et al. 2011]. Results presented in the latter studies are the best of our knowledge, but it should be highlighted that the characterization of both luminance and chrominance channels by using the 2-D auto regressive (AR) models with the same order, seems a limitation because the two channels are independent. However, it is noteworthy to precise that the results shown in Table 5 for the two studies consider only the case when authors use the spatial structure and chrominance information without adding the color information (3-D color histograms) in the characterization. We note that there is a slight performance for this later for DB5 database and for [Verdoolaege et al. 2012] for DB4 database.

We also remark from the same table the superiority of the characterization in LC color spaces compared to the RGB one except for DB4. This is theoretically confirmed since considering the image as two separable luminance and chrominance information is more intuitive and correlated with the human perception than the RGB space which considers the color as an addition of the primary colors red, green and blue. Here, we also note that [Verdoolaege et al. 2012] considered uniquely the dependency between R, G and B components (inter-channel only),

and that [Kwitt et al. 2011] considered the R, G and B dependency besides the inter-band one (for the three components). As we stressed in the introduction, inter-channel dependencies leads to high-dimensionality, due to the huge dependence between R, G and B components, and to the fact that RGB representation offers no separation between spatial structure and chrominance.

Another contribution of our approach is the definition of a closed form expression for the Rao distance in the variety of copula based probability densities [START_REF] Maliani | [END_REF]]. As it is known, deriving a similarity measure in the case of copula based multivariate models is a challenging task. For this, in [START_REF] Kwitt | [END_REF]] and [START_REF] Sakji-Nsibi | [END_REF], authors proposed a Montecarlo approximation of the KL divergence. But, this approach is computationally expensive and is not deterministic, since the KL divergence differs depending on the random number generation. An alternative approach was proposed in [Kwitt et al. 2011], using the maximum likelihood (ML) selection rule as a similarity measure. This significantly reduces the execution time in comparison with the Monte-carlo based similarity measurement. A more attractive expression of the KL divergence is provided in [START_REF] Lasmar | [END_REF], where the separability between the copula space and the marginal space is used to derive a closedform and fast KL based similarity measurement. However, KL divergence does not satisfy properties of symmetry and triangular inequality, and thus is not a distance in the right sense of the parametric manifold. Table 6 shows average percentage classification along with the execution time using both similarity measures (Rao and KL) for the textures of DB2. We recall that according to equation 37, we need expression of the Rao distance for marginals and the used copula. Expressions of the Rao distance for Weibull, Gamma and vonMises can be found respectively in [START_REF] Maliani | [END_REF]], [START_REF] Reverter | [END_REF]] and [Ceolini et al. 2010]. We can see that the proposed similarity measure slightly outperforms the ML based and Monte-Carlo based approach in term of percentage classification. However, the improvement is more significant in term of runtime which is 9 times and 12 times less than the ML based and Monte-carlo based approaches respectively. We note that the experiments were done using Matlab environment on an HP Compaq equipped with an Intel Core 2 Duo CPU at 3GHZ workstation Finally, we want to validate experimentally the hypothesis made in introduction about dimensionality reduction when we consider multi-modeling in LC color spaces. We compare our approach to the full dependence based approach of subsection 2.1. Table 7 shows performance of the classifier for multi-models Lab1 and HSV1 with full feature vectors length. It can be seen from these results

that the high parametrisation does not improve performances, but on contrary it makes them decrease by about 10%. Separability between luminance and chrominance helps then to avoid redundancy in characterization which is the case in RGB where spatial structure resides in the three color components.

Conclusion

We proposed in this letter a multi-model based approach for color textures characterization in the LC color spaces. This approach takes into account the separability between luminance and chrominance. We especially took care of the circular/linear nature of chrominance components. We also addressed the dimensionality reduction issue for the study of the full dependence in term of intra-band, inter-band and inter-channel relation between coefficients color component subbands. Results on five wide range of databases show the superiority of the proposed approach in comparison with the existing unimodeling approaches.

We know that our method as any compared methods, presents some limitations when using large databases and also against rotation and illumination effects.

Furthermore, considering the runtime issue, results show the effectiveness of the Rao geodesic distance in comparison with the Kullback-leibler divergence.

For 2 -

 2 D wavelet transforms, each scale usually has several subbands corresponding to the different orientations. Let us consider K subbands, where K = S × O with S the number of scales and O the number of orientations.

  It refers to the dependence between coefficients of subbands from different color channel, considering a fixed orientation k. The observation matrix is as follows:

2. 4

 4 .1. Multivariate spatial model for luminance (M L ) As a first step for reducing dimension of the full dataset D, the spatial structure along intra-and inter-band dependence are studied only for luminance via the model M L . For a subband coefficient at the position s, intra interband dependence is represented by the vector l s = (l 1,s , • • • , l m,s , • • • , l d,s ), where d = m×K, m is the neighborhood size and K is the number of wavelet subbands.

  Bayes classifier to assign the most probable class to a textured image. If we consider the equiprobability of each class, the Bayes classifier consists in finding the maximum of the likelihood such as: θ b is the vector of parameters associated to the class b. The vector θ b is estimated during the learning step from the set of texture samples representing the class b. Considering equation 22, the decision of the classifier is provided by the numerical evaluation of the likelihood which is time consuming because of evaluating equation 21 from the three color channel data. In order to avoid this expensive step, a direct comparison between parametric vectors, i.e. one

  T ,x b ) . The advantage of using the softmax principal is the smoothed version of the maximum, since if a b a j for all j = b, then p(C b |x T ) 1 and p(C j |x T ) 0. The proposed approach is thus based on the selection of a similarity measure, i.e L i (x T , x b ) exhibiting a closed form in terms of model parameters. Various probabilistic measures of similarity propose this kind of properties such as the Kullback-Leibler or Hellinger measures. All these possibilities have been unified

  2) d , Σ 2 ), we can compute the Rao geodesic distance as the sum of the Rao distances of the Gaussian distribution and the Rao distances between marginals [El Maliani et al. 2011]:

•

  Lab2: refers to CopMWbl for M L and CopMGam for M Cr in Lab color space. • Lab3: refers to CopMGam model for M L and the bivariate model {Gamma,Weibull} for M Cr in Lab color space. • HSV1: refers to CopMGam for M L and the bivariate model {vonMises,Weibull} for M Cr in HSV color space. (vonMises for the Hue marginal, and Weibull for the Saturation marginal) • HSV2: refers to CopMWbl for M L and the bivariate model {vonMises,Gamma} for M Cr in HSV color space. • HSV3: refers to CopMGam for M L and and the bivariate model {Weibull,Gamma} for M Cr in HSV color space.

Table 1 :

 1 Averaging measure of dependence on the Vistex database, using Pearson, Kendall, Spearman and mutual information for 40 color textures from the Vistex database in the Lab color space.

			Intra-band			Inter-band			Inter-channel
		l	cr 1	cr 2	l	cr 1	cr 2
	Pearson	0.47 0.041 0.0372 0.33 0.037 0.041	0.003	0.003	1
	Kendall	0.45 0.038 0.037 0.28 0.037 0.038	0.002	0.002	1
	Spearman 0.32 0.026 0.025 0.25 0.025 0.026	0.001	0.001	1
	MI	4.46 1.56	1.56	4.20 1.03	1.04	0.74	0.75	5.88

{l, cr 1 } {l, cr 2 } {cr 1 , cr 2 }

Table 2 :

 2 Average of parameters (considering wavelet subbands for 40 color textures from Vistex database), supposing that subbands of the three color components are described by the same marginale (Weibull or Gamma) in the case of Lab, HSv and RGB color spaces

		Weibull	Gamma
		scale	shape scale shape
	L	11.47	1.54	4.97	2.18
	a	3.21	1.63	1.26	2.41
	b	2.20	1.62	0.92	2.39
	H 26.83	1.08	27.15	1.29
	S	0.12	1.53	0.05	2.16
	V 0.0081	1.52	0.03	2.21
	R 31.51	1.53	13.83	2.16
	G 29.91	1.52	13.33	2.13
	B 31.25	1.52	13.85	2.14

Table 3 :

 3 Notations, descriptions and pdf expressions for the used models

	Model	Description	pdf
	CopMGam	Gaussian copula based multi-
		variate Gamma (all marginals
		are Gamma)

Table 4 :

 4 Performances of different versions of the multi-model in Lab and HSV color spaces. L and M Cr enable us to consider different parametric modeling for marginals between luminance and chrominance channels. As already mentioned,

		DB1			DB2		DB3	DB4		DB5
		Sn	Sp	Sn	Sp	Sn	Sp	Sn	Sp	Sn	Sn
	Lab1 97.5 0.96 98.3 0.98 89.7 0.88 58.2 0.55 77.6 0.75
	Lab2	94.5 0.92	95	0.93 84.6 0.84 54.3 0.52 74.7 0.73
	Lab3	93.3 0.92 93.8 0.90	82	0.82 53.6 0.52 73.4 0.71
	HSV1 95.6 0.94 96.7 0.94	86	0.86 55.5 0.54 76.9 0.75
	HSV2 93.1 0.93 93.4 0.93 83.1 0.83 53.6 0.53 73.9 0.73
	HSV3 90.4 0.88 90.8 0.88 80.5 0.82 51.3 0.49 70.4 0.68
	models M

Table 5 :

 5 Average percentage classification (Sn) using the best multi-models in Lab and HSV in comparison with existing approaches.

		DB1 DB2 DB3 DB4 DB5
	Lab1 (proposed method)	97.5 98.3 89.7 58.2 77.6
	[Qazi et al. 2011], in Lab	97.2 96.5	88	52.1 78.5
	HSV1 (proposed method)	95.6 97.8	86	55.5 76.9
	[Qazi et al. 2010], in IHLS	95.4 97.4 84.1 49.3 77.1
	[Verdoolaege et al. 2012], in RGB 92.2 94.7 78.5 59.6 75.1
	[Kwitt et al. 2011], in RGB	91.5 94.1 79.8 55.4 74.9

Table 6 :

 6 Average percentage classification Sn and runtime (in minutes) on DB2, using Rao distance, ML based divergence and Monte-carlo based KL divergence as similarity measure for Lab1 and HSV1 multi-models.

		Rao		KL ML		KL MC
	Sn	Runtime	Sn	Runtime	Sn	Runtime
	Lab1 98.3 4.9 min	98	40.5 min 97.2	60 min
	HSV1 96.7 6.5 min 95.9 44.3 min 96.2	63 min

Table 7 :

 7 Classifier performances for full feature vector (spatial structure for luminance and

	chrominance)
	DB1	DB2	DB3
	Sn	Sp	Sn	Sp	Sn	Sp
	Lab1 86.9 0.83	90	0.88 75.8 0.78
	HSV1 84.5 0.82 89.6 0.88 72.2 0.70

as the training set and the black half is used as the testing set. Hence, the training procedure will account for non uniformity of the original images. Every color component of each subimage was normalized by subtracting its mean and dividing by its standard deviation, and then decomposed using a 2scale dual-tree complex wavelet transform (DTCWT) [Kingsbury 1998] with a Q-shift(14,14) tap filter. In addition to its rich directional analysis, shift invariance and low redundancy properties, the DTCWT was chosen for its reduced computational time.

Quantitative evaluation of performance

Evaluation of the classification performance of each approach consists on two criteria, namely the percentage of classification and the precision. We recall that for each class of textures:


Iñaki Esnaola

Samir M Perlaza

H Vincent Poor

Oliver Kosut

DECENTRALIZED MMSE ATTACKS IN ELECTRICITY GRIDS

Keywords: Data-injection attacks, state estimation, game theory, decentralized attacks

Decentralized data-injection attack construction with minimummean-square-error state estimation is studied in a gametheoretic setting. Within this framework, the interaction between the network operator and the set of attackers, as well as the interactions among the attackers, are modeled by a game in normal form. A novel utility function that captures the trade-off between the maximum distortion that an attack can introduce and the probability of the attack being detected by the network operator is proposed. Under the assumption that the state variables can be modeled as a multivariate Gaussian random process, it is shown that the resulting game is a potential game. The cardinality of the corresponding set of Nash Equilibria (NEs) of the game is analyzed. It is shown that attackers can agree on a data-injection vector construction that achieves the best trade-off between distortion and detection probability by sharing only a limited number of bits offline. Interestingly, this vector construction is also shown to be an NE of the resulting game.

INTRODUCTION

The introduction of advanced sensing and communication infrastructure in electricity grids enables the implementation of applications and services envisioned in the smart grid paradigm but it also opens the door to cyber-security threats [START_REF] Hossain | Smart Grid Communications and Networking[END_REF]. In this paper, data-injection attacks [START_REF] Liu | False data injection attacks against state estimation in electric power grids[END_REF] against electricity grids are studied in a decentralized setting. The attack construction is formulated within a Bayesian framework in which the statistical structure of the state variables is exploited. With growing data mining and analysis capabilities provided by modern computing, it is reasonable to assume

The work of S. M. Perlaza was supported in part by the European Commission under Marie Sklodowska-Curie Individual Fellowship No. 659316 (CYBERNETS). The work of H. V. Poor was supported in part by the U.S. National Science Foundation under Grants CMMI-1435778 and ECCS-1549881.

that network operators can learn the statistical structure of the system and incorporate it into a model of the underlying stochastic process governing the network [START_REF] Ozay | Machine learning methods for attack detection in the smart grid[END_REF]. Following the approach of [START_REF] Cui | Coordinated data-injection attack and detection in the smart grid: A detailed look at enriching detection solutions[END_REF], [START_REF] Kosut | Malicious data attacks on the smart grid[END_REF] and [START_REF] Huang | Defending false data injection attack on smart grid network using adaptive CUSUM test[END_REF], the state variables are modeled as a multivariate Gaussian process whose second order moments are available to the attacker and the operator. The rationale for this is to consider a worst-case scenario for the operator. The validity of the Gaussian assumption about the distribution of the state variables is corroborated with real data for the case of low voltage distribution systems in [START_REF] Genes | Recovering missing data via matrix completion in electricity distribution systems[END_REF].

Given the complexity and extent of most electricity grids, it is plausible to think of scenarios in which several attackers intrude upon a network at different locations. In this scenario, in which multiple attackers are present and/or limited communication is available among different instantiations of the same attacker, raises the notion of distributed attacks. Distributed attack and detection strategies are investigated in [START_REF] Cui | Coordinated data-injection attack and detection in the smart grid: A detailed look at enriching detection solutions[END_REF] and [START_REF] Ozay | Sparse attack construction and state estimation in the smart grid: Centralized and distributed models[END_REF]. The decentralized system with different actors operating over a large number of processes poses a suitable framework for the exploration of game theoretic techniques. A comprehensive account of smart grid services and applications that can be tackled with game theory is given in [START_REF] Saad | Gametheoretic methods for the smart grid: An overview of microgrid systems, demand-side management, and smart grid communications[END_REF]. In [START_REF] Esnaola | Equilibria in data injection attacks[END_REF], centralized data-injection attacks are studied in a game theoretic setting in which the operator performs least squares estimation. Attack constructions that aim to manipulate market prices are modeled as a zero-sum game in [START_REF] Esmalifalak | Bad data injection attack and defense in electricity market using game theory study[END_REF]. However, the case in which several attackers disrupt the state estimation process in an uncoordinated way is still not well understood. Furthermore, the impact of making the statistical structure of the state variables available to attackers in decentralized settings has not been studied either. These issues are considered in this paper.

The next section describes the system model, including the estimation and detection procedures. The decentralized case and the properties of the resulting game are analyzed in Section 3. The paper ends with with concluding remarks in Section 4.

SYSTEM MODEL

Let x ∈ R N be a vector containing the state variables of a power system with N buses [START_REF] Abur | Power System State Estimation: Theory and Implementation[END_REF]. Assuming linearized system dynamics with M measurements corrupted by additive white Gaussian noise, the measurement vector y o ∈ R M is given by

y o = Hx + z, (1)
where H ∈ R M ×N is the Jacobian of the linearized system dynamics around a given operating point and z ∼ N (0, σ 2 I M ) is additive white noise with power spectral density σ 2 . The data-injection attack a is an M -dimensional deterministic vector introduced by an external attacker with the aim of disrupting the state estimation procedure. The attacker interferes with the measurements and modifies the observation model to

y = Hx + z + a, (2)
where y ∈ R M are the measurements available to the network operator that have been corrupted by the data-injection attack.

MMSE State Estimation

The aim of the network operator is to obtain an estimate x of the state vector x using the measurements y. In practical state estimation settings, linear estimators are preferred due to their simplicity, and thus the estimation procedure reduces to x = Ly, given a linear estimation matrix L. In the case in which the operator knows the underlying random process governing the state of the network, a common performance criterion for the estimation is to minimize the mean square error (MSE). In this case, the network operator uses a linear estimation matrix M that is the unique solution to the following optimization problem:

M ∆ = arg min L∈R M ×M E 1 N x -Ly 2 , (3)
where the expectation is taken with respect to the distributions of x and z. Under the assumption that the network state vector

x follows an M -dimensional real Gaussian distribution with zero mean and covariance matrix Σ xx , the minimum MSE (MMSE) estimation matrix is

M = Σ xx H T (HΣ xx H T + σ 2 I) -1 , (4)
and the MMSE estimate of the state vector x is

xMMSE = My. (5)
The aim of an attacker is to choose a data-injection vector a ∈ R M in order to obstruct the ability of the network operator to estimate the variables without being detected. Note that the impact of the data-injection vector a on the estimate xMMSE is quantified by the second term on the right-hand side of the following equality:

xMMSE = M(Hx + z) + Ma. (6)
The term Ma is referred to as the excess distortion induced by the attack vector a and is denoted by

x a = Ma = Σ xx H T (HΣ xx H T + σ 2 I) -1 a. (7)

Attack Detection

As a part of grid management, a network operator performs bad data detection to identify corrupted measurements. This operation can be cast as a hypothesis testing problem with hypotheses

H 0 :
There is no attack (8)

H 1 : Measurements are compromised. ( 9
)
Assuming the operator knows that x ∼ N (0, Σ xx ) it can obtain the joint probability density function of the measurements, y, and the state variables x. From ( 2) and the random model adopted to describe the state variables, it follows that the observations y are realizations of an M -dimensional real Gaussian random variable with covariance matrix

Σ yy = HΣ xx H T + σ 2 I, (10)
and mean a when there is an attack; or zero mean when there is no attack. The resulting hypothesis testing problem compares the following hypotheses:

H 0 : y ∼ N (0, Σ yy ) (11) H 1 : y ∼ N (a, Σ yy ). (12)
A worst case scenario approach is assumed for the attackers, namely, the operator knows the attack vector, a, used in the attack. However, the operator does not know a priori whether the grid is under attack or not, which accounts for the need for an attack detection strategy. That being the case, the optimal attack detection strategy for the operator is to perform a likelihood ratio test L(y, a) with respect to the measurement vector y. Under the assumption that state variables follow a multivariate Gaussian distribution, the likelihood ratio is given by

L(y, a) = f N (0,Σyy) (y) f N (a,Σyy) (y) = exp 1 2 a T Σ -1 yy a -a T Σ -1 yy y , (13)
where f N (µ,Σyy) is the probability density function of a multivariate Gaussian random variable with mean µ and covariance matrix Σ. Therefore, either hypothesis is accepted by evaluating the inequalities

L(y, a) H0 ≷ H1 τ, (14)
where τ ∈ [0, ∞) is tuned to set the trade-off between the probability of detection and the probability of false alarm. The average probability that the network operator is unable to detect the attack vector a is

P ND (a) = E 1 {L(y,a)>τ } , (15)
where the expectation is taken over the state variables x and the noise z, and 1 {•} denotes the indicator function. Note that under these assumptions, y is a Gaussian random vector with mean a and covariance matrix Σ yy . Thus, the probability P ND (a) of a vector a being a successful attack, i.e., a nondetected attack is given by [START_REF] Poor | An Introduction to Signal Detection and Estimation[END_REF] P ND (a) =

1 2 erfc   1 2 a T Σ -1 yy a + log τ 2a T Σ -1 yy a   . ( 16
)

DECENTRALIZED ATTACK CONSTRUCTION

This section describes the decentralized construction of datainjection attacks when K attackers are present. Let K = {1, . . . , K} be the set of attackers that can potentially perform a data-injection attack on the network. Let also C i be the set of sensors that attacker i controls. Assume that C 1 , . . . , C K are proper sets and form a partition of the set M of all measurement sensors. The set A k of data attack vectors a k that can be injected into the network by attacker k ∈ K is of the form

A k = {a k ∈ R M : (a k ) j = 0 for all j / ∈ C k , a T k a k ≤ E k }. (17)
The constant E k < ∞ represents the energy budget of attacker k. Let the set of all possible sums of the elements of i and A j be denoted by A i ⊕ A j . That is, for all a ∈ A i ⊕ A j , there exists a pair of vectors (a i , a j ) ∈ A i × A j such that a = a i + a j . Using this notation, let the set of all possible data-injection attacks be denoted by

A = A 1 ⊕ A 2 ⊕ . . . ⊕ A K , (18)
and the set of complementary data-injection attacks with respect to attacker k be denoted by

A -k = A 1 ⊕ . . . ⊕ A k-1 ⊕ A k+1 ⊕ . . . ⊕ A K . (19)
Given the individual data-injection vectors a i ∈ A i , with i ∈ {1, . . . , K}, the global attack vector a is

a = K i=1 a k ∈ A. (20)

Choice of Utility Function

The aim of attacker k is to corrupt the measurements obtained by the set of meters C k by injecting an error vector a k ∈ A k that maximizes the damage to the network, i.e., the excess distortion, while avoiding the detection of the global datainjection vector a. Clearly, all attackers have the same interest but they control different sets of measurements, i.e., C i = C k , for any pair (i, k) ∈ K 2 . For modeling this behavior, attackers use the utility function φ : R M → R, to determine whether a data-injection vector a k ∈ A k is more beneficial than another a k ∈ A k given the complementary attack vector

a -k = i∈{1,...,K}\{k} a i ∈ A -k (21)
adopted by all the other attackers. The utility function φ is chosen considering the fact that an attack is said to be successful if it induces a non-zero distortion and it is not detected. Alternatively, if the attack is detected no damage is introduced into the network as the operator discards the measurements and no estimation is performed. Hence, given a global attack a, the distortion induced into the measurements is 1 {L(Hx+z+a,a)>τ } x T a x a . However, attackers are not able to know the exact state of the network x and the realization of the noise z before launching the attack. Thus, it appears natural to exploit the knowledge of the first and second moments of both the state variables x and noise z and consider as a metric the expected distortion φ(a) that can be induced by the attack vector a:

φ(a) = E 1 {L(Hx+z+a,a)>τ } x T a x a , ( 22
) = P ND (a) a T Σ -1 yy HΣ 2 xx H T Σ -1 yy a, ( 23
)
where the expectation is taken over the distribution of the state variables x and the noise z. Note that under these assumptions of global knowledge, this model considers the worse case scenario for the network operator. More specifically, the result presented in this section corresponds to a conservative case in which the attackers inflict the most harm.

Game Formulation

The benefit φ(a) obtained by attacker k not only depends on its own data-injection vector a k , but also on the data-injection vectors a -k of all the other attackers. This becomes clear from the construction of the global data-injection vector a in (20), the excess distortion x a in [START_REF] Genes | Recovering missing data via matrix completion in electricity distribution systems[END_REF] and the probability of not being detected P ND (a) in [START_REF] Esnaola | Decentralized maximum distortion MMSE attacks in electricity grids[END_REF]. Therefore, the interaction of all attackers in the network can be described by a game in normal form

G = K, {A k } k∈K , φ . (24)
Each attacker is a player in the game G and it is identified by an index from the set K. The actions player k might adopt are data-injection vectors a k in the set A k in (17). Given a vector of data-injection attacks a -k , player k aims to adopt a datainjection vector a k such that the expected excess distortion

φ(a k + a -k ) is maximized. That is, a k ∈ BR k (a -k ) , (25)
where the correspondence BR k :

A -k → 2 A k is the best response correspondence, i.e., BR k (a -k ) = arg max a k ∈A k φ (a k + a -k ) . (26)
From this perspective, a game solution that is particularly relevant for this analysis is the Nash equilibrium (NE) [START_REF] Nash | Equilibrium points in n-person games[END_REF].

Definition 1 (Nash Equilibrium) The data-injection vector a is an NE of the game G if and only if it is a solution of the fix point equation

a = BR (a) , (27)
with BR : A → 2 A being the global best-response correspondence, i.e.,

BR (a) = BR 1 (a -1 ) ⊕ . . . ⊕ BR K (a -K ) . (28)
Essentially, at an NE, attackers obtain the maximum benefit given the data-injection vector adopted by all the other attackers. This implies that an NE is an operating point at which attackers achieve the highest expected distortion induced over the measurements. More importantly, any unilateral deviation from an equilibrium data-injection vector a does not lead to an improvement of the average excess distortion. Note that this formulation does not say anything about the exact distortion induced by an attack but rather it characterizes the average distortion. This is because the attack is chosen under the uncertainty of the state vector x and the noise term z.

The following proposition highlights an important property of the game G in (24).

Proposition 1 The game G in (24) is a potential game.

Proof: The proof follows from the observation that all the players have the same utility function φ [START_REF] Monderer | Potential games[END_REF]. Thus, the function φ is a potential of the game G in (24) and any maximum of the potential function is an NE of the game G.

In general, potential games [START_REF] Monderer | Potential games[END_REF] possess numerous properties that are inherited by the game G in (24). These properties are detailed by the following propositions Proposition 2 The game G possesses at least one NE.

Proof: Note that φ is continuous in A and A is a convex and closed set; therefore, there always exists a maximum of the potential function φ in A. Finally from Lemma 4.3 in [START_REF] Monderer | Potential games[END_REF], it follows that such a maximum corresponds to an NE.

Achievability of an NE

The attackers are said to play a sequential best response dynamic (BRD) if the attackers can sequentially decide their own data-injection vector a k from their sets of best responses following a round-robin (increasing) order. Denote by a (t) k ∈ A the choice of attacker k during round t ∈ N and assume that attackers are able to observe all the other attackers' datainjection vectors. Under these assumptions, the BRD can be defined as follows.

Definition 2 (Best Response Dynamics) The players of the game G are said to play best response dynamics if there exists a round-robin order of the elements of K in which at each round t ∈ N, the following holds:

a (t) k ∈ BR k a (t) 1 + . . . + a (t) k-1 + a (t-1) k+1 + . . . + a (t-1) K . ( 29
)
From the properties of potential games (Lemma 4.2 in [START_REF] Monderer | Potential games[END_REF]), the following proposition follows.

Lemma 1 (Achievability of NE attacks) Any BRD in the game G converges to a data-injection attack vector that is an NE.

The relevance of Lemma 1 is that it establishes that if attackers can communicate in at least a round-robin fashion, they are always able to attack the network with a data-injection vector that maximizes the average excess distortion.

Cardinality of the set of NEs

Let A NE be the set of all data-injection attacks that form NEs. The following theorem bounds the number of NEs in the game.

Theorem 1 The cardinality of the set A NE of NEs of the game G satisfies

2 |A NE | C • rank(H) (30)
where C < ∞ is a constant that depends on τ .

Proof: The proof of Theorem 1 can be found in [START_REF] Esnaola | Decentralized maximum distortion MMSE attacks in electricity grids[END_REF].

Theorem 1 shows that the set of attackers only need to share at most log 2 C • rank(H) bits for coordinating an attack. Note that this exchange of information needs to take place only once and can be done offline. Interestingly, the authors have not been able to obtain any example in which |A NE | > 2. In view of the numerical evidence the following conjecture is postulated. Conjecture 1 suggests that exchanging one bit of information is enough for the attackers to agree on the construction of an attack vector that leads to an NE.

CONCLUSION

This paper has established that decentralized attack construction strategies are feasible in a setting in which multiple attackers have limited communication. The trade-off between MMSE estimation and attack detection has been used to propose a novel utility function. This utility function gives rise to a game theoretic formulation that models the interaction among multiple attackers in the system. We have shown that the resulting game is a potential game, proved the existence of at least two NEs, and shown that the number of NEs is upper bounded by a finite number that depends on the network characteristics. Therefore the attackers cannot agree on an attack construction that will lead to an NE without coordination. Nonetheless, exchanging a finite number of bits offline enables the attackers to agree on a strategy that leads to an NE.

Conjecture 1

 1 The set of NEs of G satisfies |A NE | = 2.