Jose Rubio-Hernan

Luca De Cicco

SAMOVAR Joaquin Garcia-Alfaro
email: joaquin.garcia_alfaro@telecom-sudparis.com

Adaptive control-theoretic detection of integrity attacks against cyber-physical industrial systems

The use of control-theoretic solutions to detect attacks against cyber-physical industrial systems is a growing area of research. Traditional literature proposes the use of control strategies to retain, eg, satisfactory close-loop performance, as well as safety properties, when a communication network connects the distributed components of a physical system (eg, sensors, actuators, and controllers). However, the adaptation of these strategies to handle security incidents is an ongoing challenge. In this paper, we survey the advantages of a watermark-based detector against some integrity attacks as well as the weaknesses against other attacks. To cover these weaknesses, we propose a new control and security strategy that complements the watermark-based detector. We validate the detection efficiency of the new strategy via numeric simulation. Experimental results are also presented by using a laboratory testbed based on supervisory control and data acquisition industrial protocols.

INTRODUCTION

In an effort to reducing complexity and costs, industrial infrastructures are being upgraded with novel computing, communication, and interconnection capabilities. The goal is to build and operate them as easily and reliably as possible. Current trends, such as smart factories and intelligent supply chains mainly rely on sensors, actuators, and controllers that are networked to communicate with each other by using cyber-physical industrial protocols. [START_REF] Hespanha | A survey of recent results in networked control systems[END_REF][START_REF] Kim | Cyber-physical systems: a perspective at the centennial[END_REF] However, this opens challenging security threats. [START_REF] Wu | A survey on the security of cyber-physical systems[END_REF] The recently coined cyber-physical security term refers to mechanisms that address security issues associated to these environments. The use of inadequate cyber-physical security mechanisms can have an adverse effect in smart industries. [START_REF] Falliere | W32 stuxnet dossier[END_REF][START_REF] Weyer | Towards industry 4.0 -Standardization as the crucial challenge for highly modular, multi-vendor production systems[END_REF][START_REF] Lee | A cyber-physical systems architecture for industry 4.0-based manufacturing systems[END_REF] These issues place the study of cyber-physical security mechanisms as a hot research topic.

Given the control-theoretic nature of cyber-physical industrial systems, the control community is actively working to adapt traditional control strategies to detect faults and errors, toward detectors of malicious attacks. [START_REF] Salt | Sistemas de control basados en red modelado y diseño de estructuras de control[END_REF][START_REF] Heemels | Periodic event-triggered control for linear systems[END_REF][START_REF] Han | Stochastic event-triggered sensor schedule for remote state estimation[END_REF] Motivated by the same objectives, we present in this paper a solution that combines 2 different control strategies to handle integrity attacks against cyber-physical systems. The new solution combines an event-triggered control strategy together with a challenge-response watermark-based detector. We present the theoretical approach underlying our proposed strategy as well as its algorithmic implementation. We validate our proposed approach via numerical simulations, complemented by experimental results obtained using a laboratory testbed based on supervisory control and data acquisition (SCADA) industrial protocols.

This paper is organized as follows. Section 3 provides introductory background. Section 4 presents our detection strategy. Section 5 validates the approach via numerical simulations. Section 6 complements the simulations with experimental results obtained with a laboratory testbed based on SCADA protocols. Section 2 surveys related work. Section 7 concludes this paper.

RELATED WORK

Security of cyber-physical systems is drawing a great deal of attention recently. [START_REF] Corman | NITRD cyber-physical security panel[END_REF] Solutions focused on control approaches for the detection of cyber-physical attacks are the research axis more closely related to this paper. This axis is the one that explicitly considers the interconnection between cyber and physical control domains in networked control systems. Recently, the control system community started to study the security of cyber-physical systems both under the methodological point of view and from a more technological standpoint by looking at particular problems arising in, eg, smart grids. Concerning the methodological aspects, several studies have proposed to adapt classical frameworks to handle security issues in networked control systems.

In a cyber-physical replay attack, the adversary aims at carrying out malicious actions without holding a priori knowledge about system model dynamics. Most of the other cyber-physical attacks reported in the literature require that the adversary gathers some system knowledge. For example, to execute a dynamic false-data injection attack, [START_REF] Mo | False data injection attacks against state estimation in wireless sensor networks[END_REF] the adversary has to have a perfect knowledge of the plant's behavior. Likewise, to execute a covert attack, [START_REF] Smith | Covert misappropriation of networked control systems: presenting a feedback structure[END_REF] it is necessary to have some knowledge about the model of the plant and the controller. The adversaries handled in this paper fit those aforementioned models defined in the works of Garone et al, [START_REF] Mo | False data injection attacks against state estimation in wireless sensor networks[END_REF] and Smith, [START_REF] Smith | Covert misappropriation of networked control systems: presenting a feedback structure[END_REF] ie, it is assumed that they aim at obtaining knowledge about the model underlying the behavior of the plant and the controllers prior attacking the system.

Concerning detection mechanisms, previous work in the related literature considers the adaptation of fault detection systems to detect cyber-physical attacks. [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF][START_REF] Teixeira | A secure control framework for resource-limited adversaries[END_REF][START_REF] Pasqualetti | Cyber-physical security via geometric control: Distributed monitoring and malicious attacks[END_REF] In particular, Mo et al [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] and Mo and Sinopoli [START_REF] Mo | Secure control against replay attacks[END_REF] proposed the use of the watermark-based detection by adapting traditional failure detection mechanisms to detect replay attacks. This watermark is added to the control measurements to verify using the detection mechanism that the sensor measurements are not replayed measurements, ie, the control measurements with the watermark have to be correlated with the sensor measurements. Miao et al [START_REF] Miao | Stochastic game approach for replay attack detection[END_REF] improved the performance of this detection mechanism. They present a suboptimal algorithm against replay attack, using a stochastic game approach. This suboptimal algorithm combines the watermark-based detector proposed by Mo et al and the zero-sum stochastic game proposed by Zhu and Başar. [START_REF] Zhu | Dynamic policy-based IDS configuration[END_REF] In the same way, Do et al [START_REF] Do | A statistical method for detecting cyber/physical attacks on SCADA systems[END_REF] formulated the attack detection problem as a transient changes detection problem in stochastic-dynamical systems. This detector is based on the knowledge of the behavior of the system and its stochastic variations to detect data manipulation. Using this detection control algorithm allows to protect the system against attacks able to perform malicious actions in a short period of time, of which goal is to disrupt safety-critical applications. An alternative approach is taken by Pasqualetti et al, [START_REF] Pasqualetti | Cyber-physical security via geometric control: Distributed monitoring and malicious attacks[END_REF] where the authors propose to use the theory of geometric control to model cyber-physical systems attacks. In line with the aforementioned approach, this paper focuses on the interconnection between control strategies and a watermark detector to handle the integrity attacks.

BACKGROUND

Cyber-physical attacks

Some of the abbreviations used hereinafter are summarized in Table 1. The coined cyber-physical industrial system term integrates a physical infrastructure and a cyber framework in an effort to automate physical processes, as well as reducing complexity and costs of traditional industrial systems. These traditional systems are composed of sensors, actuators, and other devices that interact with the physical processes. The technology evolution brings these systems toward a combination between a physical layer, which encompasses the physical framework, and a cyber layer that encompasses the communication and computation framework, [START_REF] Genge | A system dynamics approach for assessing the impact of cyber attacks on critical infrastructures[END_REF] eg, by using SCADA technologies. Cyber-physical systems combine control strategies through the communication network that are used to manage the physical layer with some cyber strategies and used to monitor the physical layer. The cyber layer has to follow some constraints driven by control strategies that also appear as constraints at the physical layer. Attacks targeting both layers are often referred to as cyber-physical attacks. In a nutshell, cyber-physical attacks aim at affecting physical processes of a networked control system by manipulating network information at the cyber layer.

Teixeira et al [START_REF] Teixeira | A secure control framework for resource-limited adversaries[END_REF] proposed a taxonomy of cyber-physical attacks based on the resources of the adversaries. For instance, the knowledge of the adversary about the system is an important resource to assess the nature of the attacks. Based on the degree of the adversary knowledge, the attacks may succeed at violating some system properties, such as availability and integrity, as well as to obtaining operational information about the system itself to cover the attacks from detection. Mo et al [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] addressed the traditional replay attacks, in which the adversary does not hold any a priori knowledge about the system dynamics. They also address injection attacks, as a variant of the previous type, in which the adversary injects false data or alterations from legitimate data.

Those aforementioned attacks are also addressed by Teixeira et al [START_REF] Teixeira | A secure control framework for resource-limited adversaries[END_REF] but classified in a more granular manner (depending on the techniques to inject the malicious data) as surge attacks, bias injection attacks, and geometric attacks. In order to perform such attacks, the adversary must hold some partial knowledge about the physical processes underlying the target system. Smith [START_REF] Smith | Covert misappropriation of networked control systems: presenting a feedback structure[END_REF] presented a more powerful attack type, referred to as covert attacks, in which the adversary must perfectly know the system dynamics, to evade detection. To acquire the system knowledge, adversaries are assumed to learn and identify as much information as possible from the network or physical data. Rubio-Hernan et al [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF][START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF] differentiated between parametric and nonparametric adversaries, depending on the identification techniques used to obtain such system knowledge.

In terms of attack detection, several techniques have been proposed in the literature. Wang et al [START_REF] Wang | SRID: State relation based intrusion detection for false data injection attacks in SCADA[END_REF] proposed the use of relation graphs to identify false data from injection attacks. Arvani and Rao [START_REF] Arvani | Detection and protection against intrusions on smart grid systems[END_REF] proposed a signal-based intrusion detection model to identify anomalies in the information reported by the sensors of cyber-physical systems. Do et al [START_REF] Do | A statistical method for detecting cyber/physical attacks on SCADA systems[END_REF] presented the use of statistical detection methods, to identify data manipulation over stochastic processes. The use of watermark-based techniques are proposed by Mo et al [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] and revisited by Rubio-Hernan et al. [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF][START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF] Next, we reexamine the watermark-based techniques presented in the works of Mo et al [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] and Rubio-Hernan et al [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF][START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF] and show how to combine such techniques together with control strategy policies. We show that the resulting combination succeeds at handling the most powerful types of cyber-physical adversaries reported in the literature, ie, cyber-physical adversaries holding complete knowledge about the system dynamics.

Control strategies

Control theory is a well-known topic, where the evolution of the technology has been the main motivation to create new control policies to manage these systems, keeping the control features. Among these new technologies, we can mention the networked control systems, where the loop between the different components of the system is closed through the network. The evolution of networked control systems to deal with the industry 4.0 is defined as cyber-physical industrial system. These new defined systems need the collaboration among a very wide number of disciplines to solve the challenges in term of autonomy, reliability, usability, functionality, and cyber security. [START_REF] Baheti | Cyber-physical systems[END_REF] Figure 1 shows a simple cyber-physical industrial system diagram, where f s and f c are the frequencies used by the sensors and controller, respectively, to process data. f ca is the sampling frequency used in the communication channels between the controller and the actuators. f sc is the sampling frequency used in the communication channels between the sensors and the controller.

Among control strategies used in cyber-physical industrial systems, we focus on strategies depending on the transmission policy necessary to secure the system. We assume policies not only to identify control faults but also intentional attacks carried out by external entities. Such policies can be classified as either sampled-data control policies or event-triggered control policies.

In the case of sampled-data control policies, the system can be driven using either monofrequency sampling (ie, same sampling frequency for all the channels) or multifrequency sampling (ie, different sampling frequencies depending on the channel, such as sensor-to-controller channel, or controller-to-actuator channel [START_REF] Salt | Sistemas de control basados en red modelado y diseño de estructuras de control[END_REF] ). Figure 1, case 1, represents a monofrequency sampling scenario. Figure 1, case 2, represents a multifrequency sampling scenario.

Security in cyber-physical systems shall include the management of control properties and policies through the network, to avoid that an external entity, eg, an adversary, succeeds at harming the system. In the case of event-triggered control policies, we can differentiate between continuous event-triggered control [START_REF] Heemels | Model-based periodic event-triggered control for linear systems[END_REF] policies, where all the data received by the sensors are sent to the controller; periodic event-triggered control [START_REF] Heemels | Periodic event-triggered control for linear systems[END_REF] policies, where sensors send update data periodically; and stochastic event-triggered schedule [START_REF] Han | Stochastic event-triggered sensor schedule for remote state estimation[END_REF] policies, where sensors send update data only when there are variations in data obtained. Figure 1, cases 3, 4, and 5, represent each of these different types of policies.

Watermark-based attack detection

We address in this section the watermark-based detector proposed by Mo et al. [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] The goal of such a detector is to identify cyber-physical replay attacks against networked control industrial systems. To analyze the watermark-based detector, Mo et al used an industrial control system modeled mathematically as a discrete linear time-invariant system. This mathematical model is used to describe the dynamic behavior of the system. The system can be represented as follows:

x t+1 = Ax t + Bu t + w t ( 1
)
y t = Cx t + v t , (2)
where x t ∈ R n is the state's vector, u t ∈ R p is the control signal, y t ∈ R m is the system output, and w t ∈ R n and v t are the process noise and the measurement noise, respectively. The noises are assumed to be a zero mean Gaussian white noise with covariance Q, ie, w t ∼ N(0, Q) and R, ie, v t ∼ N(0, R), respectively. Moreover, A ∈ R n×n , B ∈ R n×p , and C ∈ R m×n are the state matrix, the input matrix, and the output matrix, respectively. Let us now define the well-known linear quadratic Gaussian (LQG) approach used as a control technique in the work of Mo et al. [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] This technique has 2 independent components. 1. a Kalman filter producing an optimal state estimation xt of the state x:

xt|t-1 = Ax t-1 + Bu t-1 xt = xt|t-1 + K t (y t -Cx t|t-1 ), (3)
where K t denotes the Kalman gain, and xt|t-1 is the a priori system state estimation. 2. a linear quadratic regulator providing the control law u t .

u t = Lx t , ( 4
)
where L denotes the feedback gain of a linear-quadratic regulator.

After describing the model of the plant, hereinafter, we present the detection scheme proposed in the work of Mo et al [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] against replay attacks. The idea is to superpose a watermark signal Δu t ∈ R p to the optimal control law u ⋄ t . The new control input u t is given by

u t = u ⋄ t + Δu t . ( 5
)
Note that the watermark signal is independent from the process noise w t and the output noise v t . To detect the adversaries, the watermark-based detector employs a well-known χ 2 detector. [START_REF] Brumback | A chi-square test for fault-detection in Kalman filters[END_REF] Figure 2 shows the overall control system equipped with the attack detector proposed in the work of Mo et al. [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] The alarm signal g t generated by the detector is defined as follows:

g t = t ∑ i=t-w+1 (r i ) T  -1 (r i ), (6)
where w is the size of the detection window,  is the covariance of input signals from the sensors, and r t = y t -Cx t|t-1 is the residue generated from the estimator at each tth time step.

LQR Controller

Detector FIGURE 2 Watermark protection in cyber-physical systems. [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF]

LQR, linear quadratic regulator

To verify if the system is under attack, g t is compared with threshold γ. If g t is equal or greater than the threshold, g t ≥ γ, then the detector generates an alarm. The aforementioned threshold is chosen with respect to the system properties, eg, during the system identification phase, to keep a proper ratio of false positives and detector efficiency.

Watermark-based attack detection against cyber-physical adversaries

The stationary challenge-response watermark-based detector proposed in the work of Mo et al [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] has been analyzed and revisited in the work of Rubio-Hernan et al. [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF][START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF] In the work of Rubio-Hernan, [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF] a new adversary, named nonparametric cyber-physical adversary is defined to show some weaknesses of using a stationary watermark in the detector. The adversary is able to use nonparametric identification techniques to identify the behavior of the system from the input and output data of the system, ie, from data sent to the actuator and data received from the sensors. This adversary bases the attack on the stationarity of the data. These stationary data allows the adversary to identify the correlation between the different events executed in the system. To handle such a vulnerability, a new watermark-based detector is proposed in the work of Rubio-Hernan et al. [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF] The new approach uses different distributions to move toward an adaptive, and nonstationary, watermark-based detector with 2 main configurable parameters: the number of distributions N and the switching frequency f = 1∕T. The new watermark is defined as follows:

Δu t = Δu (s(t,T)) t , ( 7
)
where

s ∶ N × R → ,  = {0, 1, … , N -1}
, is a static function that maps the time sample t and the switching period T to an element of the index set , defined as follows:

s(t, T) = ⌊ 1 T mod (t, NT) ⌋ , ( 8
)
where mod(x, y) is the modulo operator, and ⌊•⌋ is the floor function.

The new watermark-based detector has been analyzed in the work of Rubio-Hernan et al [START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF] against a more powerful adversary. The new adversary, as in the case of the nonparametric adversary, uses parametric identification techniques to obtain a precise version of the system model. Unlike the nonparametric adversary, the techniques used by the new adversary aim at computing some precise system parameters. This allows to obtain the behavior of the system regardless of using nonstationary watermarks. To do so, the new adversary needs to guess 2 main configuration parameters: the precise number of samples to identify the system and the order of the system, ie, the dimension of the correlation matrix between the input and the output. This new adversary is able to evade detection if those 2 aforementioned parameters are properly chosen. In the following section, we present a new strategy to handle the adversary presented in the work of Rubio-Hernan et al. [START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF]

Performance loss

The quadratic cost J of a LQG controller, with respect to the state x t and the control input u t , can be obtained using the following expression: where Γ and Ω represent positive definite cost matrices. [START_REF] Franklin | Digital control of dynamic systems[END_REF] The strategy proposed in this paper uses 2 detectors previously proposed in the works of Mo and Sinopoli [START_REF] Mo | Secure control against replay attacks[END_REF] and Rubio-Hernan et al. [START_REF] Rubio-Hernan | Revisiting a watermark-based detection scheme to handle cyber-physical attacks[END_REF] Both detectors use the theoretical methodology proposed by Mo et al, [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] to compute optimal detection, in terms of performance loss. Our approach does not increase the performance loss with respect to previous work.

J = lim n→∞ E [ 1 n n-1 ∑ i=0 ( x T i Γx i + u T i Ωu i ) ] , (9)

PIETC-WD WATERMARK-BASED DETECTION STRATEGY

In the previous section, we have seen that the watermark-based schemes are able to handle attacks carried out by adversaries with limited knowledge about the system dynamics, for instance, the ones defined as either cyber adversaries or nonparametric cyber-physical adversaries (cf Section 3.4). Nevertheless, it fails at detecting those adversaries with enough knowledge about the system dynamics, defined as parametric cyber-physical adversaries (cf Section 3.4). In this section, we present a new detector scheme, hereinafter denoted as periodic and intermittent event-triggered control watermark detector (PIETC-WD). This new detector aims at detecting cyber adversaries as replay attack and the cyber-physical adversaries presented in previous section.

Our scheme consists of a local controller located in each sensor and a remote controller creating a distributed controller (cf Figure 3). The cooperation between the local and the remote controller allows us to create an intrusion detection policy to capture integrity attacks. The local controllers manage the dynamics of the plant, and the remote controller manages the closed-loop system to ensure the system against integrity attacks. Notice that our new scheme requires an additional controller, together with the sensors, which must have enough computation power to process data estimations, eg, to predict errors between environmental and estimated data. The actuators do not require additional computational power. Nevertheless, during the time between 2 consecutive events, they must keep the last data received from the remote controller.

To perform with our scheme, it is necessary to define communication policies among the sensors, the actuators, and the remote controller. We define 2 communication policies for ensuring the system: (1) periodic communication policy, which the communication from the sensors to the remote controller is periodical, with a T sc = 1∕f sc period, and from the remote controller to the actuators, with a T ca = 1∕f ca period; and (2) intermittent communication policy, which allows for sending data from the sensors to the remote controller if the local controller produces an alarm. Notice that T sc cannot be equal to T ca to avoid that an intermittent communication takes place while the periodic communication is being sent. Definition 1. Periodic and intermittent event-triggered control watermark detector is a detector strategy with distributed control tasks. On one hand, the sensors control the system periodically, using their local controllers and a local watermark-based detector. [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] On the other hand, the remote controller uses the estimation error received from each sensor to periodically generate the control inputs. The remote controller also controls the closed-loop communication with an intermittent watermark.

The PIETC-WD strategy uses essential properties from cyber-physical systems, to detect the attacks. The 3 main properties used in our work are as follows: (1) stability of closed-loop systems; (2) controllability of the system, to get full reachability; and (3) observability, to correlate events from the system. The periodic communication policy of PIETC-WD uses these 3 properties to verify if the system works properly, with regard to the information (eg, measurements) provided by the sensors. The intermittent communication policy verifies the stability of the closed loop from the controller standpoint, to avoid manipulation from external agents.

The communication algorithms, used to perform the PIETC-WD strategy, are shown in Algorithms 1 and 2. Algorithm 1 shows the local controller implementation, placed in the sensors, whose input is the data obtained from the physical system, y i t , and its output is the residue of the local controllers, r i t (with a challenge response watermark), if the alarm is disabled; or the value obtained by the sensors, y i t , if the alarm is enabled. Algorithm 2 shows the remote controller implementation, whose input is the data sent by the sensors, data sc ; and its output is the control inputs sent by the controller, data ca , and the alarm value, alarm c . More information about the controllers and the communication policies is provided in the sequel. The following notation will be used hereinafter:

U = {u 0 , u T ca , … , u 𝜅T ca } control input variables Y = {y 0 , y T sc , … , y ηT sc } observed measurements,
where αT ca = T sc with α, η and 𝜅 ∈ N.

Local controller design

The local controller is located in the sensors and uses a watermark to verify that the dynamics of the system is correct. Each sensor has a local controller with an LQG approach (cf Section 3.3). We denote the local controller in each sensor by i ∈ {0, 1, … , N -1}, where N is the number of sensors in the system. This controller adds a watermark to the sensor measurement before sending the residue to the remote controller

y (i) t = y ⋄(i) t + Δy (i) t ( 10
)
r (i) t = y (i) t -C i x(i) t|t-1 , ( 11
)
where y ⋄(i) t is the sensor measurement, Δy (i) t is the watermark added by the local controllers, and r (i) t is the residue sent to the remote controller to compute the control input u t = [u (0) , … , u (N-1) ]. Notice that the new sensor measurement y (i) t is computed after verifying that y ⋄(i) t is the correct sensor measurement.

Remote controller design

The remote controller receives periodically the residue of each sensor, r (i) t , and computes these residues using the LQG approach (cf Section 3.3) to obtain the state estimation

xt = xt|t-ϱ + K t (r t ) (12)
and

r t = r τ ∀t ∈ [τ, τ + (α -1)T ca ] , ( 13
)
where r t is a vector generated by all the residues of the sensors and r τ with τ ∈ [t ∶ mod(t, T sc ) = 0] is the periodical residues' vector. In Equation 12, xt|t-ϱ is defined as

xt|t-ϱ = A ϱ xt|t-ϱ + ϱ ∑ j=1 A j-1 Bu ϱ-j+1 , (14)
where ϱ ∈ {1, … , α}.

Equation 12 works properly if the data used by each sensor are independent among them or the system has only one sensor. Otherwise, if the data received by sensors are correlated, it is necessary to add a rectification for each sensor's residue to consider the cross-correlation among the watermarks added by local controllers. We can compute this rectification for t ∈ [T 0 , T 0 + T -1] as follows:

r rect t (i) = r t (i) -C i Δx i t|t-T , ( 15
)
where r rect t (i) is the rectified value respect to different sensors. In addition, Δx i t|t-T is the rectification.

Δx i t|t-T = t-T 0 -1 ∑ j=0 ( (A i + B i L i ) j+1 K i (ϵ r ) ) , ( 16
)
where ϵ r = r t-1-jr (i) t-1-j . We can define the control inputs vector, u t , as follows:

u t = L ( xt|t-ϱ + K t r rect t ) = L ( xt|t-ϱ + K t (r * t + Δy t ) ) , ( 17
)
where r * t is the vector of the residues before adding the watermark and their respective rectification, and Δy t is the vector generated by all the sensors' watermarks.

The watermark used intermittently by the remote controller is added to the control inputs. The controller adds a watermark with probability β. Denoting λ t = 1 or 0 as indication function whether the watermark is added or not, we assume that λ ′ s are independent and identically distributed (iid). Bernoulli random variables with E[λ t ] = β.

The intermittence of the watermark communication allows us to define the watermark behavior as a nonstationary distribution. This watermark, Δu t (cf Equation 5), permits us to detect if the closed loop is being manipulated. It is worth noting that Δu t is a stochastic signal with the same variance as Δy t .

Periodic communication policy

The periodic communication policy is managed by the sensors. The sensors add the watermark in the measurements received by the plant and send residue r t to the remote controller. The remote controller uses these residues to generate the control inputs sent to the actuators. The actions of these actuators produce change in the state of the plant that is captured by the sensors. If the real states differ from the state estimated by the sensors, then the sensors will switch from periodic communication policy to intermittent communication policy (cf Section 4.3).

In order to validate the proposal, let us assume that an attack is started at time T 0 , and we compute residue

r (i) for t ∈ [T 0 , T 0 + T -1] r (i) t = y ′(i) t -C i x(i) t|t-T , ( 18
)
where y ′(i) t is the sensor measurement sent to the controller by the adversary. Moreover, it is easy to show that the following holds:

x(i) t|t-T = x′(i) t|t-T +  t-T 0 i ( x(i) T 0 |T 0 -1 -x′(i) T 0 |T 0 -1 ) + t-T 0 -1 ∑ j=0 (  j (A i + B i L i )K i (ϵ Δy ) )
where ϵ Δy = (Δy (i) t-1-j -Δy ′(i) t-1-j ), x′(i) is the local estimated state for each sensor when the system is under attack, and [START_REF] Mo | Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs[END_REF] Substitution of 19 in 18 yields

 i = (A i + B i L i )(I i -K i C i ) is a stable matrix.
r (i) t = y ′(i) t -C i x′(i) t|t-T ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟ First term -C i  t-T 0 i ( x(i) T 0 |T 0 -1 -x′(i) T 0 |T 0 -1 ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ Second term -C i t-T 0 -1 ∑ j=0 (  j (A i + B i L i )K i (ϵ Δy ) ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ Third term .
Let us consider separately the 3 terms in the equation written above: the first term follows the same distribution of (y t -C i x(i) t|t-1 ); since  i is asymptotically stable, ie, all its eigenvalues are inside the open unit disk of the complex plane, the second term converges exponentially to zero. In fact, the entries of  t-T 0 i converge exponentially fast to zero. The third term, under attack, is not equal to zero, since Δy (i) t ≠ Δy ′(i) t , and the adversary is detected; for a cyber adversary viewpoint, the measurements of the sensors change all the time and replay measurements are not accepted; likewise, a cyber-physical adversary is not able to obtain the system model using the methodology proposed in Section 3.4. The parametric cyber-physical adversary model, H at , using the autoregressive with exogenous input identification approach, [START_REF] Natke | System identification: Torsten Söderström and Petre Stoica[END_REF] is computed as follows:

H at = f at (u t , r t , y t , v t ), (19)
where f at is a linear function that map the input u t , as a function of t to the output (residue r t , or measurement y t ) using a process noise, v t , to measure a stochastic error.

Assuming that the real model is, H = f(u t , y t , v t ), a linear function. Furthermore, the model shown by the system is H = f s (r t , y t , v t , u t , Δu t , Δy t ), a non linear function, because we add an Δy t in the measurements of the sensors periodically and an Δu t in control input, as stochastic event, avoiding to follow a linear function to describe the system. For this reason, a powerful cyber-physical adversary is detected.

Intermittent communication policy

The aforementioned periodic communication policy is managed by the sensors. The sensors produce an alarm if g t ≥ γ. When a sensor produces an alarm, this information is sent immediately to the remote controller. The affected sensor sends the real measurement to the remote controller to perform a second verification. This alarm is activated in a sensor if the control input has been manipulated by an external entity, a problem occurs in the system or the remote controller adds the watermark in the control input.

As we have defined before, a remote controller can generate an intermittent communication to verify that the closed loop is not broken. For a SISO (single-input, single-output) system or SIMO (single-input, single-output) system, the remote controller follows the next steps.

• The controller uses a low-pass filter after computing the control input, u t , to remove the variations, and generate an optimal control input, u * t .

u * t = L ( xt|t-ϱ + K t r * t ) (20)
• The remote controller adds the owns watermark.

u t = L ( xt|t-ϱ + K t r * t ) + Δu t ( 21
)
The local controllers located in the sensors receive the modification. They generate an alarm and send the measurement of the sensors to the remote controller. Then, the remote controller verifies if the measurement of the sensors is correct using the χ 2 detector. If the detector value is under the threshold, the system works correctly. Otherwise, if this value is over the threshold, the remote controller generates an alarm. In both cases, the next control input is computed as follows:

u t+T ca = L ( xt|t-ϱ + K t r t ) , ( 22
)
where residue r t follows Equation 13. For a MISO (multiple-input, single-output) system or a MIMO (multiple-input, multiple-output) system the remote controller follows the next steps:

• The controller computes the control input:

u t = L(x t|t-ϱ + K t r t ) = [ u (1) , … , u (i) , … , u (n) ] , ( 23
)
where n is the number of control inputs. • Then the new control input vector is generated as follows:

u t { u t = [ u (1) , … , u (i) , … , u (n) ] ∀i ≠ j u (i) t = u * (i) t + Δu (i) t i = j ,
where u * ( j) t

is the control input of the actuator, j, filtered to remove the variations. After removing the variation, the remote controller adds its watermark. It is worth noting that j is the actuator chosen by the controller to send the watermark.

As we have shown before, if the detector value, generated from the real data sent by the sensor, is under the threshold, the next control input is computed as follows:

u t+T ca = L(x t|t-ϱ + K t r t ). (24)
Otherwise, the remote controller generates an alarm. Briefly, when the remote controller receives a measurement from a sensor, if a watermark Δu has not been sent, then the remote controller creates an intrusion alarm. Otherwise, if a watermark has been added to the control input, the controller verifies if this alarm is produced by the watermark. If the residue generated between the real measurements of the sensors and the estimation is under the threshold, the remote controller sends the control input generated before adding the watermark. However, if the residue is over the threshold, it means that an external entity is into the closed loop, and an alarm is activated.

In order to validate our claims, let us assume the following attack in the communication channel between the sensor and the controller after the controller sends a control input with a watermark. It is started at time T 0 , and the remote controller includes the remote watermark, Δu t at time

T 1 ∈ [T 0 , T 0 + T -1]. We compute the residues r t for t ∈ [T 0 , T 0 + T -1]: r t { r t = [ r (1) , … , r (i) , … , r (n) ] t ∈ [T 0 , T 1 ] r t = y ′ t -Cx t|t-T t ∈ [T 1 , T 0 + T -1] . ( 25
)
On one hand, it is worth noting that the attack is not detected between T 0 and T 1 , since the adversary is able to replay or insert a correct residue without being detected. On the other hand, for t ∈ [T 1 , T 0 + T -1], it is easy to show that the following holds:

xt|t-T = x′ t|t-T +  t-T 1 ( xT 1 |T 1 -1 -x′ T 1 |T 1 -1 ) + t-T 1 -1 ∑ j=0 (  j B ( Δu t-1-j -Δu ′ t-1-j )) . ( 26
)
Substitution of ( 26) in ( 25) yields:

r t = y ′ t -Cx ′ t|t-T ⏟⏞⏞⏞⏟⏞⏞⏞⏟ First term -C t-T 1 ( xT 1 |T 1 -1 -x′ T 1 |T 1 -1 ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ Second term -C t-T 1 -1 ∑ j=0 (  j B ( Δu t-1-j -Δu ′ t-1-j
))

⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟

Third term .

The first term follows the same distribution of (y t -Cx t|t-1 ); the second term converges exponentially to zero. Since the third term is not equal to zero, Δu t ≠ Δu ′ t , the adversary is detected; from the cyber adversary viewpoint, the measurements of the sensors change all the time and replay measurements are not accepted; likewise, the cyber-physical adversary is not able to obtain the system model.

New parametric cyber-physical adversary

In this section, we present a new parametric cyber-physical adversary with the knowledge about the new detector strategy, to evaluate new detection strategy. This attacker has knowledge about the new communication policies and the existence of the local and the remote watermarks. Nevertheless, the new adversary does not know the watermark covariances, the controller's parameters used to obtain the correct error between data, and neither the moment when the remote controller forces an intermittent communication.

The new adversary could be able to detect the correlation model between the inputs and the outputs of the plant. This adversary can force the sensors' intermittent communication with malfunction control inputs, and mislead the controller with replay error data to obtain the model. Nevertheless, this adversary is not able to know when the communication is periodic or intermittent, since the attacker does not know when the remote control sends the watermark added to the control inputs, which generate the intermittent communication. The intermittent communication does not change the communication between the remote controller and the actuators, but produces an intermittent communication between the sensors and the remote controller, necessary to verify the closed loop.

Briefly, the new adversary is able to attack the integrity of the system. Nevertheless using the PIETC-WD strategy, the adversary is detected by the controllers of the sensors. The remote controller detects the attack when it verifies the behavior of the closed loop. The adversary cannot avoid the alarm in the sensors (local controller). Nevertheless, the attacker can cut off the communication between the sensors and the remote control misleading the remote controller with correct residues (eg, replay residues). Moreover, to avoid the alarm in the remote controller, the adversary can switch between sending the measurements of the sensors or the residues, but the adversary has a great probability to be detected. We validate the PIETC-WD strategy against the new parametric cyber-physical adversary in the next section.

NUMERICAL VALIDATION

This section presents a numerical validation of the strategy. We have simulated a multiple-input-multiple-output system, which represents a cyber-physical system with 4 inputs and 4 outputs. To define mathematically this system using the model shown FIGURE 4 The dynamics of the states vector in the plant under a parametric attack. The attack starts at t = 700s in Section 3.3, we use the following matrices * : We present in Figure 4 the dynamic of the plant and how the adversary can modify this dynamic to disrupt the system. Figure 5 shows how the detectors present in the sensors are able to detect attacks using the periodical communication policy and a watermark managed by sensors. Figure 5A,B shows the dynamic of the states estimated by the local controllers placed in sensors and the remote controller placed in control center, respectively. Figure 5C,D shows the dynamic of the alarm in sensors (local detectors) and in the control center (remote detector). We show in these Figures an attack, which is detected by the local detectors, but whose alarm's dynamic is not transmitted to the control center. The sensors' alarms are stopped by a cyber-physical adversary that takes the control of the system and sends a dynamic alarm correlated with the data sent from the remote controller (control center) to the actuator. This action allows the adversary attack the system without be detected.

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -0.
Figure 6 presents the control/security strategy using the periodical and intermittent communication policies. As shown in Figure 6, it is visible that using both policies allows the system to detect the attack. Figure 6C shows the local detector that is able to detect a disruption in the dynamic of the plant, thanks to the periodical communication policy. In Figure 6D, we can see that the remote detector (placed in the control center) is able to detect that the closed loop is broken, using the remote controller's watermark and the intermittent communication policy.

Our strategy handles integrity attacks at the physical level, by using the watermark sent by the local controller, ie, placed in the sensors (cf Figure 6A), as well as integrity attacks carried out through the network, via the closed-loop verification conducted by the remote controller, using stochastic metrics to trigger the watermark (cf Figure 6B). The states estimated by the local and remote controllers (cf Figure 6A,B) are handled and mitigated by the system, reporting the events to security operators (eg, in terms of alerts) and activating the appropriate corrective mechanisms.

*The covariance matrices are equal to Q = 0.2I and R = I.

EXPERIMENTAL RESULTS

We present some experimental results obtained with a laboratory testbed. The testbed implements the detection strategy and the adversary models presented in this paper. The implementation is conducted using Modbus 30 and DNP3, [START_REF] Curtis | A DNP3 protocol primer. A basic technical overview of the protocol[END_REF] 2 well-known SCADA protocols used in cyber-physical industrial systems.

SCADA Testbed

SCADA is a general term that encompasses well-defined types of field devices, such as: (1) master terminal units and human machine interfaces, located at the topmost layer and managing all communications; (2) remote terminal units (RTUs) and programmable logic controllers (PLCs), controlling and acquiring data from remote equipment and connecting with the master stations; and (3) sensors and actuators. The design of our testbed includes all those aforementioned field devices: controllers, sensors, actuators, PLCs, and RTUs. The devices are distributed across several nodes in a shared network combining Modbus 30 and DNP3. [START_REF] Curtis | A DNP3 protocol primer. A basic technical overview of the protocol[END_REF] One or various field devices can be embedded into a single node. From a software standpoint, the controller never connects directly to the sensors. Instead, it is integrated in the testbed as a PLC node, with eventual connections to some other intermediary nodes. To run our experiments, we created a distributed network where the control center communicates with the physical devices through the RTUs. The communication between the control center and the RTUs uses the protocol DNP3 and the communication between the RTUs and the physical devices uses Modbus. The testbed can be expanded to handle additional SCADA protocols, as well as to connect complementary devices, such as additional PLCs and RTUs.

From a data transmission standpoint, the devices connected to the testbed use different sampling frequencies. The testbed is able to handle several PLCs. To avoid overloading 1 channel with all the possible registers of the PLCs, different ports are designated to isolate the communication between separated PLCs. DNP3 commands perform an integrity scan that gathers all the data from the PLCs in case several PLCs were being handled in the same channel. This avoids variables of the PLCs causing communication overheads. (C) Control center detector results not using intermittent policy.

(D) Control center detector results using intermittent policy. A sample cyber-physical industrial scenario implemented in our testbed is shown in Figure 7. In the Figure, the solid and pattern gray squares represent 2 mobile agents that move from 2 spatial coordinates (represented by the the solid and pattern gray circles), and vice versa. Each mobile agent controls the dynamics of its movements. The 2 agents are coordinated by a SCADA controller, to avoid spatial collisions. Further information, as well as some live demonstration videos, is available at http://j.mp/TSPScada.

Experimentation and results

We present the results of applying the detection strategy, assuming attackers under the parametric adversary model. The attacker tries to evade the detection techniques reported in this paper. In a first scenario, the system uses only a periodic communication policy using the watermark detector placed at the sensors (cf Figure 8A,C). When the attacker successfully infers the system parameters, the attack is not reported by the remote controller. Nevertheless, local alarms placed at the sensors report the attack. In a second scenario, the system uses also the intermittent policy. In such a case, the attack is detected (cf Figure 8B,D) by the remote controller. In Figure 8B, we may observe some peaks (cf, t = [2.1s; 14.4s; 15.3s; 17.6s; 20.4s; 21.0s; 22.3s; 26.2s; 37.8s; 56.8s; 66.8s]), representing the local detector's reaction in the presence of a watermark sent by the remote controller -used to verify that the closed loop works properly, as defined in Section 4.3.

To compare the watermark-based detector mechanism and our detection strategy, several repetitions of the experiment were orchestrated. Using the results obtained from the testbed experiments, we report in Table 3 the percentage of false positives, false negatives, and detection ratio. We also show the average time that the detector placed at the control center, takes to launch the alerts. Notice that the workload in the sensors, due to adding the local controller in our detection strategy, does not significantly increase the performance loss of the system. The increase, in terms of memory and processing time, is shown in Table 2.

Regarding the results shown in Table 3, we can emphasize that: (1) a system that uses only the watermark-based detector mechanism against a parametric attacker has a lower detection ratio, of about 12%. That is possible since the attacker can evade the detection process if they succeed at properly identifying the system attributes. This result is expected, as suggested by the theoretical and simulation-based conclusions presented in the work of Rubio-Hernan [START_REF] Rubio-Hernan | Event-triggered watermarking control to handle cyber-physical integrity attacks[END_REF] ; and (2) a system that uses the strategy proposed in this paper has a higher detection ratio, of about 65.00%. In this scenario, the detection ratio increases, thanks to the distributed detector proposed in Section 4, confirming the theoretical and simulation-base result reported in the previous sections. Moreover, the false negative ratio significantly decreases, from 88.60% to 35.00%. In terms of false positives, both scenarios show similar data, but the PIETC-WD strategy generates about 2.0% more. The time between the beginning of the attack and the moment, when the attack is detected by the remote controller, also increases, since the detection watermark managed by the remote controller follows a stochastic law.

CONCLUSION

We have addressed security issues in cyber-physical industrial systems. We have presented an adaptive strategy that detects adversaries that are able to acquire knowledge about the system dynamics, prior to starting their attacks to get control over the inputs and measurements of the system. In addition to industrial systems, any other cyber-physical scenario, modeled in terms of sensors, actuators and feedback control, can benefit from the approach. Computer networks complemented by physical systems, such as Internet of things scenarios, may also benefit from the approach. We have validated the construction using numerical validation and a laboratory testbed based on SCADA industrial protocols. We have shown that the strategy is able to detect the attacks with high detection ratio.

FIGURE 3

 3 FIGURE 3 Cyber-physical system diagram with the new control security strategy

(

  A) Estimated states in the local controller. (B) Estimated states in the remote controller. (C) Sensor detector results. (D) Control center detector results.

FIGURE 5

 5 FIGURE 5Numeric simulation results, using only periodical policy. Attacks start at t = 700s. A,B, The dynamics of the states vector estimated in the local and remote controller, respectively, under the same scenarios. C,D, The dynamics of the alarm signal g t produced in the local and remote detector, respectively, under the same scenarios

(

  A) Estimated states in the local controller. (B) Estimated states in the remote controller. (C) Sensor detector results. (D) Control center detector results.

FIGURE 6

 6 FIGURE 6Numeric simulation results, using periodical and intermittent policies. Attacks start at t = 700s. A,B, The dynamics of the states vector estimated in the local and remote controller, respectively, under the same scenarios. C,D, The dynamics of the alarm signal g t produced in the local and remote detector, respectively, under the same scenarios

FIGURE 8

 8 FIGURE 8 Experimental testbed results. The horizontal solid line represents the threshold. The vertical dotted line represents the moment when the attack Peaks on the left side of the vertical dotted line represent false positives. A,B, Local detector values g t (placed within the sensors) without and with intermittent policies (ie, remote watermarks) under a parametric attack. C,D, Remote detection values (placed at the control center) without and with intermittent policies under a parametric attack

TABLE 1

 1 Abbreviations used in this paper

	Acronyms	Descriptions
	CPS	Cyber-physical system
	SCADA	Supervisory control and data acquisition
	NCS	Networked control system
	IT	Information technology
	ETC	Event-triggered control
	PETC	Periodical event-triggered control
	LTI	Linear time invariant
	LQG	Linear quadratic Gaussian

PIETC-WD

Periodical and intermittent event-triggered control and watermark detector

PLC

Programmable logic controller

RTU

Remote terminal units

ARX

Autoregressive with exogenous input

  Cyber-physical system diagram. Case 1: f ca = f sc . Case 2: f ca ≠ f sc . Case 3: f ca = f c , f sc = f s , and f s = f c . Case 4: f ca =

	Plant	Sensors
	Network
	Controller
	FIGURE 1 f c p a	, f sc =	f s p s	,
	such that p a , p s ∈ N. Case 5: f sc and f ca →stochastic events

TABLE 2

 2 Memory and processing time increase, due to the PIETC-WD strategy implementation

	Time	RAM Memory
	Sensor's workload increase 17.13 ms	5.23%
	Abbreviation: PIETC-WD, periodic and intermittent event-triggered
	control watermark detector.

TABLE 3

 3 Detector performance results

		Using Only Watermark Using PIETC-WD
		Detector Mechanism	Strategy
	Detection ratio	12.00%	65.00%
	Average detection time	6.08 s	7.30 s
	False negatives	88.60%	35.00%
	False positives	1.35%	3.5%
	Abbreviation: PIETC-WD, periodic and intermittent event-triggered control watermark
	detector.

ACKNOWLEDGEMENTS

The authors acknowledge support from the Cyber CNI Chair of the Institut Mines-Télécom (cf http://chaire-cyber-cni.fr/). The chair is supported by Airbus Defence and Space, Amossys, EDF, Orange, La Poste, Nokia, Société Générale and the Regional Council of Brittany. It has been acknowledged by the Center of excellence in Cybersecurity.


Bruno Borloz
email: borloz@univ-tln.fr

-Bernard Xerri
email: xerri@univ-tln.fr

Subspace SNR maximization: the constrained stochastic matched filter

Keywords: detection, subspace method, reduced-rank method, signal-to-noise ratio maximization, matched filter, matched subspace

In this paper, we propose a novel approach to perform detection of stochastic signals embedded in an additive random noise. Both signal and noise are considered to be realizations of zero mean random processes whose only secondorder statistics are known (their covariance matrices).

The method proposed, called Constrained Stochastic Matched Filter (CSMF), is an extension of the Stochastic Matched Filter itself derived from the Matched Filter. The CSMF is optimal in the sense that it maximizes the Signalto-Noise Ratio in a subspace whose dimension is fixed a priori.

In this paper, after giving the reasons of our approach, we show that there is neither obvious nor analytic solution to the problem expressed. Then an algorithm, which is proved to converge, is proposed to obtain the optimal solution.

The evaluation of the performance is completed through estimation of Receiver Operating Characteristic curves. Experiments on real signals show the improvement brought by this method and thus its significance.

I. Introduction

This paper deals with the problem of detecting a stochastic signal (like a transient signal for example) embedded in an additive random noise.

Throughout this paper, all the signals will be real and discrete (time samples, pixels of images, ...) and represented with vectors of R N .

The method proposed here consists in a linear filtering called (for reasons explained later) "Constrained Stochastic Matched Filter" (CSMF). This method gives, for an integer value p (1 ≤ p < N ), among all the p-dimension subspaces, the one where the Signal-to-Noise Ratio (SNR) is maximum: the CSMF is optimal for this criterion. This is a reduced-rank method (a projection) under constraint (the constraint being the a priori knowledge of the dimension p) [START_REF] Borloz | Estimation, détection et classification par maximisation du rapport signal-à-bruit : le filtre adapté stochastique sous contrainte[END_REF].

The SNR is invariant in a p-dimension subspace: it does not depend on the basis chosen to describe the subspace. An important consequence of this invariance of the SNR w.r.t. the basis is that the simplest basis, say an orthonormal one, can usefully be chosen: moreover in such a basis the mathematical expression of the SNR is simple to obtain and will simplify later calculations.

In this paper we show that there is neither immediate nor obvious way to find the optimal p-dimension subspace: then we propose an algorithm and its the convergence to the good solution is proved.

The performances of the method and the comparisons with other methods are performed through the Receiver-Operating-Characteristic (ROC) curves giving the Probability of Detection P D w.r.t. the Probability of False Alarm P F A . Nevertheless, this paper gives no demonstration that ROC curves are better for a predicted value of p: we only observe, with results obtained from numerical simulations, that there exists a value of p for which the ROC curve is the best one.

Let us also note that our model is not a parametric one. The only knowledge is the covariance matrices of the random signals .

A. Problem Formulation

Let us consider an observation x ∈ R N . Two hypotheses can be formally stated (detection problem): this measurement was produced by ambient noise n alone or by a signal s embedded in this noise, respectively:

H 0 : x = n H 1 : x = s + n
The objective is to decide between these hypotheses. Our model will not be a parametric one.

The assumptions of our model are the following : 1) s and n are realizations of zero mean random processes.

2) The covariance matrices of s and n, respectively A and B, are supposed to be known, full rank and different.

3) s and n are uncorrelated, not necessarily Gaussian, and their Probability Density Functions (PDF) are unknown.

Two kinds of error are possible: the missing of the signal and the false alarm. A trade-off (highlighted by the ROC curves) must be found between a small average number of misses and a small average number of false alarms.

When the PDF of the signals are known, the key quantity to compute is the Likelihood Ratio (LR) L(x) which must be compared to a threshold determined according to a criterion such as the minimization of the probability of error, the maximization of P D when P F A is fixed a priori (Neyman-Pearson criterion) [START_REF] Van | Trees, Detection, Estimation and Modulation Theory[END_REF]- [START_REF] Hero | Signal Detection and Classification[END_REF].

When the PDF are unknown, L(x) cannot be calculated. This is why we take into consideration methods based on SNR maximization. Furthermore the CSMF method takes place among numerous currently known reduced-rank techniques which have been proposed (Section II is a survey of some SNR maximization and reduced-rank methods justifying the approach of our method).

B. Why the CSMF ?

When the PDFs are known (e.g. Gaussian), the Likelihood Ratio Test (LRT) is an optimal test which leads to compare a value to a threshold. For Gaussian signals, the test can easily be written as a sum of N terms (then onto the whole space of the used signals) depending on the observation vector x:

log L(x) = Λ(x) = N i=1 λ i 1 + λ i (u t i x) 2 . ( 1
)
where the λ i (λ 1 ≥ λ 2 ≥ ... ≥ λ N ) and u i are respectively the eigenvalues and corresponding eigenvectors of B -1 A (see Section II-C). In fact, the eigen-elements of B -1 A naturally appear when trying to maximize the output SNR of a linear filter h: this output SNR can be written ρ = h t Ah h t Bh . The maximal value of ρ, noted ρ max , is obtained for h = u 1 the eigenvector associated with the largest eigenvalue λ 1 of B -1 A: this filtering consists in a projection of the signal x onto E u1 , and then it is easy to verify that ρ max = λ 1 . This method is called "Stochastic Matched Filter" (SMF) [START_REF] Cavassilas | Extension de la notion de filtre adapté. Contribution à la détection de signaux courts en présence de termes perturbateurs[END_REF].

When signals are not Gaussian, we can continue to use Λ(x) which is no longer the log of the LR. This expression has no reason to be optimal, and experimental results show, first, that a truncation of this sum to p terms can improve the ROC curves and next that there exists an optimal value of p for which the ROC curve is the best one. This truncation (p < N ) is expressed as follows

p i=1 λ i 1 + λ i (u t i x) 2 ,
and can be seen as a projection of x onto a p-dimension subspace E † p : E † p is spanned by {u i } i=1,...,p . This method, called "Extended Stochastic Matched Filter" (ESMF) could be wrongly interpreted as a SNR maximization method: in fact it does not maximize the SNR in a p-dimension subspace but a weighted sum of output SNRs, each of them after a projection onto E ui for i = 1, ..., p (see Section II-C).

The method proposed in this paper is naturally inferred from these remarks concerning the output SNR maximization and the projection onto a subspace of dimension two or higher; therefore its aim is to maximize the SNR in an aptly chosen subspace with an a priori given dimension p. The choice of p, and then the dimension of the optimal subspace searched for, is a constraint: this is why the name of "Constrained Stochastic Matched Filter" (CSMF) was given to this optimal filter. We will clearly see that the CSMF is not a simple extension of the ESMF and that the CSMF can no more be inferred from the ESMF.

However, when p = 1, the CSMF and the SMF are identical. But when p > 1, it is proved in this paper (cf. Section III-F) that the optimal space E * p cannot be simply deduced from the knowledge of either E † p or E * p-1 . Hence, it is necessary to propose an algorithm that finds the solution: this algorithm is given and is proved to converge to the solution (cf. Section IV).

Organization of the paper

In Section I we formulate the mathematical model, present the basic assumptions. Section II describes existing methods and introduces those proposed in the paper. The method is detailed in Section III and useful properties are highlighted. Section IV is dedicated to the practical determination of this subspace: an algorithm is proposed to find the optimal subspace E * p and the proof of its convergence is given. Then experimental results are presented in Section V.

In this paper, we apply the method to detection, but it could be obviously used for compression, filtering or estimation problems.

II. Overview of some existing methods

The model is those given in Section I-A.

A. The Karhunen-Loève Transform

The Karhunen-Loève transform (KLT) is a Principal Component Analysis used to tackle this model [START_REF] Cavassilas | Extension de la notion de filtre adapté. Contribution à la détection de signaux courts en présence de termes perturbateurs[END_REF]- [START_REF] Yang | RLS-based adaptive algorithms for generalized eigendecomposition[END_REF] when noise is white (B = σ 2 n I) or absent; it provides the best approximation, in the sense that it minimizes a mean square error (MSE) for a stochastic signal under the condition that its rank is fixed and is used for example for data compression or filtering. When noise is white, it determines the p-dimension subspace where the SNR is maximum.

But it does not consider colored noise, and therefore is not optimum even when it is used with a noise suppression filter such as the Wiener filter (which is not a reduced-rank method). The SMF, a Generalized Eigen Decomposition (GED) introduced by Cavassilas-Xerri [START_REF] Cavassilas | Extension de la notion de filtre adapté. Contribution à la détection de signaux courts en présence de termes perturbateurs[END_REF] will be detailed in Section II-C: it performs a two-stage operation (prewhitening and KLT), but is shown to be not optimal in terms of maximization of the SNR. GED is a major problem in many modern information processing applications (adaptive filtering, blind source separation [START_REF] Pokharel | Principe Recursive complex BSS via generalized eigendecomposition and application in image rejection for BPSK[END_REF], ...) and fast algorithms to estimate and track the principal generalized eigenvectors have been developed [START_REF] Yang | RLS-based adaptive algorithms for generalized eigendecomposition[END_REF].

We will show that the CSMF can be seen like an extension of the KLT and the SMF for the problems we are interested in. Some authors want to find an optimal linear data compression method in the presence of noise: the Proportional KLT (PKLT) applies an oblique projection operator onto a subspace S (dim(S) = p) along a subspace L (S and L are both unknown) [START_REF] Yamashita | Relative Karhunen-Loeve transform[END_REF]. This operator P naturally maximizes a ratio of powers. Solving this problem without any constraint concerning the rank of P naturally leads to an impossibility.

The first part of their work shows how justifiable it is to take an interest in the maximization of the SNR in a subspace.

The CSMF proposed in this paper solves this problem by adding a constraint, the rank of the subspace to project data onto.

B. Parametric Models

For detection with reduced-rank methods, many authors have worked on a parametric model of the following form:

x = Hθ + n = p i=1 θ i h i + n (H is a N × p matrix)
. The useful signal s = Hθ is a stochastic signal constrained to lie in the signal subspace, the p-dimension subspace spanned by the known modes h i , with mode weights or gains, the entries θ i of θ. This model in an extension of those used for the Matched Filter (MF) detector, matched to a signal that is assumed to lie in a 1-dimensional subspace (i.e. H = h 1 is the deterministic signal to detect).

It is noteworthy that our model is very different, because the covariance matrix A is full rank and then the signal subspace has dimension N .

When the noise is Gaussian, the output of the MF provides a sufficient statistic for any LRT for detection. The knowledge of s and the second-order statistics of n is necessary to derive the corresponding MF. For p > 1 the MF detector is no more efficient and is extended to the Matched Subspace Detector (MSD) [START_REF] Scharf | Matched Subspace Detectors[END_REF]- [START_REF] Scharf | Matched and Adaptive Subspace Detectors When Interference Dominates Noise[END_REF] assuming prior knowledge of B. The MF is also named coherent MSD.

When the gains are unknown, the Generalized Likelihood Ratio Test [START_REF] Hero | Signal Detection and Classification[END_REF] takes the form of a ratio of two quadratic forms of prewhitened observations using orthogonal projections onto suitable subspaces. The statistic obtained has natural invariances (the energy of the subspace signal and the SNR are unchanged).

When B is unknown, it is obtained from training data (Adaptive Subspace Detectors) [13][14].

Numerous papers deal with the MF detector and its extensions: several reasons may imply that signal or/and noise are not exactly known (channel nonlinearities, timing jitter, non-stationarities, model uncertainties ...) [18][19].

Another problem studied (e.g. in digital communications [START_REF] Proakis | Digital Communications[END_REF]) is those of detecting a transmitted signal when one out of several known signals is transmitted. When the additive noise is white and Gaussian, the optimal detector consists of a bank of MF followed by a detector which chooses as the detected signal the one with the maximal output value. Improvements have been observed in many cases [START_REF] Eldar | Orthogonal Matched Filter Detection[END_REF].

C. The Stochastic Matched Filter and the Extended SMF

C.1 Introduction

The SMF was first introduced to detect a random signal not supposed to lie in a known subspace: furthermore, the second-order statistics of both s and n are supposed to be known [START_REF] Cavassilas | Extension de la notion de filtre adapté. Contribution à la détection de signaux courts en présence de termes perturbateurs[END_REF]. It can be seen like an extension of the MF (it provides an optimal filter since it maximizes the SNR), but also of the KLT. This problem is a generalized eigenvalue problem using the covariance matrices A and B; the filtering is a projection onto the optimal subspace spanned by the eigenvector of B -1 A with maximum eigenvalue which is also the value of the maximum output SNR. The output SNR ρ of a linear filter h can be written like a Rayleigh quotient: ρ = h t Ah h t Bh (if A and B have unit trace, ρ can be interpreted like a gain on the SNR). This problem is equivalent to solve the following generalized eigenvalue problem: Ah = λBh. The maximal value ρ max is obtained for h = u 1 the eigenvector associated with the largest eigenvalue λ 1 of B -1 A; then ρ max = λ 1 > 1: this filtering performs a projection of the signal onto E * 1 = E u1 . If {λ i } and {u i } are the eigenvalues and eigenvectors of B -1 A with λ 1 ≥ ... ≥ λ N , we easily prove that: λ i ≥ 0 can be interpreted like a gain on SNR after projection onto E ui .

{u i } is a non-orthogonal basis of R N perform- ing simultaneous diagonalization of A and B. If U = [u 1 ...u N ],
an appropriate normalization of the u i gives, with ∆ positive diagonal matrix [START_REF] Horn | Matrix analysis[END_REF]:

U t AU = ∆ U t BU = I
The interpretation of the λ i naturally leads to take into account directions of projection that could statistically contribute to a better detection, that means growing the dimension of the subspace to project data onto. Actually, it has been shown [START_REF] Cavassilas | Extension de la notion de filtre adapté. Contribution à la détection de signaux courts en présence de termes perturbateurs[END_REF] that a few other eigenvectors, the dominant one, can statistically contribute to improve ROC curves. To take a decision, we have to propose a function of x and the u i .

When the signals are Gaussian, the calculation of the logarithm of the LR leads easily to equation [START_REF] Borloz | Estimation, détection et classification par maximisation du rapport signal-à-bruit : le filtre adapté stochastique sous contrainte[END_REF]. This expression has no reason to be optimal when the signals are not Gaussian, but, according to the remarks above, this summation is shortened to p terms corresponding to significant eigenvalues and the function chosen is then:

Λ(x; p) = p i=1 λ i 1 + λ i (u t i x) 2 , ( 2
)
which is a weighted sum of (u t i x) 2 the power of the observation x after projection onto each direction u i , each weight being linked to the SNR on this direction. This method, called "Extended SMF" (ESMF), does not maximize the SNR in a p-dimension subspace, but a weighted sum of output SNRs, each of them after a projection onto E u i for i = 1, ..., p. We will denote E † p the subspace spanned by {u i } i=1,...,p .

C.2 Illustration: a practical example

To illustrate the interest of taking into account a subspace of dimension higher than one, let us apply this method on a narrow-band signal embedded in simulated underwater acoustics: the central frequency is f 0 = 3131 Hz and the spectral bandwidth of noise is B = 1260 Hz. The sampling frequency is here F e = 15423 Hz. Experiments are performed on N r = 1000 realizations of signal denoted s i (N = 21, see appendix A) and the initial SNR is about -14 dB. Envelop detection techniques could be used but, in practice, they give lower quality results.

A is calculated as follows

: A = N r i=1 s i s t i .
For these simulations, the areas of presence or not of a signal are obviously known. A detection (which can be a false alarm) is decided each time there are at least 5 consecutive points of the result function done by equation ( 2) above the threshold.

ROC curves are shown in figure 1, first for p = 1 (E † 1 ) and for the optimal value of p, say 3 (E † 3 ): The improvement brought by projecting data onto E † 3 is obvious. This example shows how the make of decision can be greatly improved by taking into consideration more than one eigenvector. Even if the SNR in E † 3 is smaller than in E † 1 , we observe that increasing the projection subspace dimension p brings an improvement that is not counterbalanced by the the decrease of the SNR.

We can also see that there is a worsening of the ROC curve for p = 4. The projection onto E † 4 will give statistically worse results than those onto E † 3 . This result confirms that there is an optimal value of p for the chosen criterion.

C.3 Conclusion

The SMF can be proved to be a two-stage operation: data whitening and then maximisation of SNR in a pdimensional subspace E † p (KLT). But, as whitening is not an optimal operation in terms of SNR, the global operation has no reason to maximize the SNR in E † p . Hence, methods which try to maximize a SNR while performing a reducedrank operation are natural when PDF are unknown. Thus, it seems natural and legitimate to ask oneself if there exists a p-dimensional subspace (p > 1) in which the SNR is maximal, with the hope that ROC curves would be improved again, and then if it is possible to find it.

III. The Constrained Stochastic Matched Filter

The method proposed in this paper has been called Constrained SMF (CSMF) because it can be seen like an extension of the SMF, naturally inferred from the remarks in previous section concerning a projection onto a subspace of dimension two or higher where the SNR is maximum.

A. Preliminary remarks and notations

A random signal s, such as in our model, is a vector of R N and can always be expressed as follows:

s = N i=1 α i x i = Xa where a = [α 1 ...α N ]
t is a vector of random variables and

X = [x 1 ...x N ] a basis of R N .
In this paper we are interested only in powers in subspaces. Let us denote A = E (ss t ) the covariance matrix and P s the power of s. Of course, P s depends only on the subspace and not on the basis used to describe it, and with no loss of generality, it is possible to consider only orthonormal bases to describe any subspace: hence X t X = I.

It readily follows that:

P s = N i=1 x t i Ax i = tr X t AX = tr AXX t = trA. (3)
Moreover, as P s = trA, we can consider, without loss of generality, only covariance matrices of trace 1.

B. SNR in a p-dimension subspace E p

Considering an integer p chosen a priori in [1, ..., N -1], let us denote E p = E {x 1 ,...,x p } the p-dimension subspace spanned by the p orthonormal vectors x 1 to x p . We will also denote E ⊥ p = E {x p+1 ,...,x N } . Then, with X p = [x 1 , ..., x p ], the projection of s onto E p along E ⊥ p has power P p :

P p = tr X t p AX p = p i=1 x t i Ax i ,
and the expression of the SNR in E p takes the form:

ρ = p i=1 x t i Ax i p i=1 x t i Bx i = tr X t p AX p tr X t p BX p . ( 4
)
The objective is to find the unknowns {x i } in order to maximize this term. The optimal subspace will be denoted E * p and the corresponding SNR ρ * . trA trB is the input SNR and if covariance matrices A and B have unit trace, ρ is in fact a gain on the SNR (and no longer the output SNR) which can be proved to be necessarily lower than the largest eigenvalue of B -1 A, say λ 1 [START_REF] Borloz | Estimation, détection et classification par maximisation du rapport signal-à-bruit : le filtre adapté stochastique sous contrainte[END_REF].

We see here the important difference with the SMF which maximizes the following expression:

ρ smf = p i=1 v t i Av i v t i Bv i ,
where the v i does not form an orthonormal basis.

Throughout the following section, we will focus our attention on equation ( 4) and try to find E * p .

C. Characterization of the optimal subspace E * p Let us consider a p-dimension subspace E p spanned by a set of p orthonormal vectors X = [x 1 ...x p ]. The expression of the SNR ρ in E p is given by the equation ( 4).

The constraints can be expressed by the p 2 relationships "x t i x j = δ ij ". Clearly, p is given and the unknowns of our problem are ρ * and the p vectors x i which must be calculated so as to maximize ρ. We are faced with an optimization problem with constraints which is usually solved using a Lagrange multipliers method. Let us define the following function:

L (X, Ω) = ρ + p i=1 p j=1 ω ij (x t i x j -δ ij ). (5)
This equation can be written:

L (X, Ω) = tr X t AX tr X t BX + tr Ω X t X -I , ( 6
)
where

Ω ≡ [ω ij ] is a p × p symmetric matrix.
A necessary condition for this value to be maximum is ∂L ∂X = 0, which means that for ρ = ρ * :

(A -ρ * B) X tr X t BX + XΩ = 0.
As B is positive definite, tr X t BX > 0 and this equation becomes

(A -ρ * B) X = XΩ 0 (7)
where Ω 0 is a p×p real symmetric matrix but not diagonal. But we can find a real orthogonal matrix Π and a real diagonal matrix ∆ * µ ≡ [µ * i ] such that Ω 0 = Π∆ * µ Π t . Then equation [START_REF] Pokharel | Principe Recursive complex BSS via generalized eigendecomposition and application in image rejection for BPSK[END_REF] becomes:

(A -ρ * B) XΠ = XΠ∆ * µ . ( 8
)
As Π is invertible, XΠ and X span the same subspace E * p . Furthermore, as Π is a real orthogonal matrix, the set of orthonormal vectors X is changed to another set of orthonormal vectors XΠ. Noting XΠ = T * = [t * 1 ...t * p ], equation ( 8) can be written

(A -ρ * B) T * = T * ∆ * µ (9)
which is an eigenvalue problem. Note that for any value of ρ, (A -ρB) is always real symmetric and then diagonalizable through a N × N unitary real eigenvector matrix. That means that E * p is spanned by a set of p orthonormal vectors which are eigenvectors of (A -ρ * B). Nevertheless, equation ( 9) is not simple to solve because, if the t * i and µ * i are naturally unknown, ρ * is unknown too. For any value ρ > 0 let us denote:

(A -ρB) t i = µ i t i , i = 1, ...N. ( 10
)
(A -ρB) is always real symmetric and the {t i } naturally form an orthonormal basis. We note that the eigenvalues µ i depend on ρ: it is easy to show, with trivial examples, that they are non nonlinear w.r.t. ρ. A simple illustration is given on figure 2 where we can see the evolution of the eigenvalues µ i w.r.t. ρ for the covariance matrices A and B given in example 1 (cf. Section III-F). D. Property of the optimal subspace E * p Equation [START_REF] Scharf | Matched Subspace Detectors[END_REF] implies that for any i, j ∈ {1, ...N }, t t j (A -ρB) t i = µ i t t j t i = µ i δ ij . Then for any subset I of {1, 2, ..., N } with cardinality p,

i∈I t t i (A -ρB) t i = i∈I t t i At i -ρ i∈I t t i Bt i = i∈I µ i , so that i∈I t t i At i i∈I t t i Bt i -ρ = i∈I µ i i∈I t t i Bt i . ( 11
)
For the optimal subspace E * p , i.e. ρ = ρ * , the left expression is null, implying that

i∈I * µ * i = 0. ( 12
)
This property will be used to find the solution in the two trivial cases 'p = 1' and 'p = N -1' (cf. Section III-E), but also to prove the convergence of the algorithm (cf. Section IV-B).

We have denoted I = I * ; it is easy to show that I * = {1, 2, ..., p} if the eigenvalues are sorted in decreasing order, say µ * 1 > µ * 2 > ... > µ * N ).

P E. Particular cases

In some particular cases, it is possible to reach the solution easily, without any sophisticated algorithm.

p = 1

The eigenvalue µ * 1 to take into account is null. Hence,

(A -ρ * B) t * 1 = µ * 1 t * 1 = 0, i.e. At * 1 = ρ * Bt * 1 .
ρ * is the largest eigenvalue of B -1 A and t * 1 its associated eigenvector. Naturally we find in this case the SMF.

p = N -1 As N -1 i=1 µ * i = 0, µ N = 1 -ρ * . Hence, (A -ρ * B) t N = (1 -ρ * ) t N , or (A -I N ) t N = ρ * (B -I N ) t N ; ρ * is the largest eigenvalue of (B -I N )
-1 (A -I N ) and t N its associated eigenvector. The N -1 vectors spanning E * N -1 are the other eigenvectors.

B = I In this case, (A -ρI) t i = µ i t i , i.e. At i = (ρ + µ i )t i . For ρ = ρ * , p i=1 µ i = 0 and then p i=1 (ρ + µ i ) = pρ = p i=1 λ A i (where λ A i is the i-th eigenvalue of A), that means ρ = 1 p p i=1 λ A i .
Hence the Karhunen-Loeve filter is a particular case of the CSMF.

F. Examples

For the following examples, we state p = 2 and search the optimal subspace E * 2 . As p = N -1, we can use the results of the previous paragraph: the SNR ρ * obtained in E * 2 is easy to calculate.

Example 1 p . This will be performed in Section IV.

A = 1 3   1.

G. Conclusion

In this section, the problem has been presented and equations have been deduced that must be solved to find the optimal p-dimension subspace. We have seen that there exists neither analytic nor obvious solution and that an algorithm must be proposed. This is the purpose of the next section.

The CSMF consists in finding the p-dimension subspace which maximizes the SNR after only a projection.

P

IV. Algorithm to find the optimal subspace E * p In light of the examples of previous section, it is required to find an algorithm to determine the optimal subspace E * p spanned by p vectors t i verifying equation ( 10) and maximizing ρ defined by

ρ = i∈I t t i At i i∈I t t i Bt i , ( 13
)
where I is a subset of p different numbers out of {1, ...N }.

It seems natural that such an algorithm should be iterative and use, at each step, these equations alternately.

A. Presentation of the algorithm

ρ * being unknown, it is reasonable to choose for the initial value of ρ, say ρ (0) , the largest eigenvalue of B -1 A. At each step n ≥ 0, we obtain the symmetric matrix A -ρ (n) B and calculate its N eigenvectors t

(n) i asso- ciated to eigenvalues µ (n)
i . Then we must choose among them the p ones t

(n) i , i ∈ I (n) /card(I (n) ) = p for which ρ (n+1) = i∈I (n) t (n) i t At (n) i i∈I (n) t (n) i t Bt (n) i (14) is maximum. These p vectors span a subspace E (n) p .
Then it is easy to calculate A -ρ (n+1) B , I (n+1) and the new subspace E (n+1) p . This process can be iterated until |ρ (n+1) -ρ (n) | < ε (see table 1). Of course, we have to prove that this algorithm converges to the good solution ρ * .

B. Study of the convergence

At step n, from ( 11) and ( 14), the variation of ρ is

ρ (n+1) -ρ (n) = i∈I (n) µ (n) i i∈I (n) t (n) i t Bt (n) i . ( 15
) Of course, if i∈I (n) µ (n) i = 0, then ρ (n+1) = ρ (n)
. Hence, as it has been proved that for ρ * , there exists a subset I * of cardinal p such as ( 12) is verified, say i∈I * µ * i = 0 , and ρ * is clearly a fixed-point of the algorithm.

Let us denote:

i∈I (n) µ (n) i = f ρ (n) . ( 16
)
Then f (ρ * ) = 0. As for any value of ρ, A -ρB is symmetrical, the {t i } span an orthonormal basis: then t t i t i = 1, ∀i. This expression can be differentiated w.r.t. ρ, leading to:

t t i ∂t i ∂ρ = 0, ∀i. ( 17
)
Description of the algorithm

ρ (0) = λ 1 the largest eigenvalue of B -1 A n = 0 1) Calculate M (n) = A -ρ (n) B 2) Compute the eigenvalues µ (n) i
and corresponding eigenvectors t

(n) i (i = 1 to N ) of M (n)
(it is possible to sort the eigen-elements so that µ

(n) 1 ≥ µ (n) 2 ≥ ... ≥ µ (n) N )
3) Find the C p N combinations of p elements out of {1, 2, ..., N }: they will be denoted

I (n) k with card I (n) k = p and 1 ≤ k ≤ C p N 4) For k = 1 to C p N Calculate ρ (n+1) k = P i∈I (n) k t (n) i t At (n) i P i∈I (n) k t (n) i t Bt (n) i (k th combination) 5) Find the maximal value of ρ (n+1) k k=1,...C p N : this maximum is denoted ρ (n+1) 6) If |ρ (n+1) -ρ (n) | < ε stop iterations, t * i = t (n) i , µ * i = µ (n) i , ρ * = ρ (n+1) Else n ← n + 1, Goto 1
-Bt i + (A -ρB) ∂t i ∂ρ = ∂µ i ∂ρ t i + µ i ∂t i ∂ρ .
From equation ( 17) and multiplying by t t i , it comes (using equation t t i (A -ρB) = µ i t t i ) :

-t t i Bt i + µ i t t i ∂t i ∂ρ = ∂µ i ∂ρ , or ∀i, ∂µ i ∂ρ = -t t i Bt i < 0. ( 18
)
Using equations ( 16) and [START_REF] Verdu | Signal Selection for Robust Matched Filtering[END_REF], equation [START_REF] Scharf | Geometries invariances, and SNR in-terpretations of matched and adaptive subspace detectors[END_REF] becomes:

ρ (n+1) = ρ (n) + i∈I (n) µ (n) i - i∈I (n) ∂µ (n) i ∂ρ = ρ (n) - f (ρ (n) ) f (ρ (n) ) , (19)
Obviously f (ρ * ) is not null: in fact, from [START_REF] Verdu | Signal Selection for Robust Matched Filtering[END_REF] we have

f (ρ) = i∈I ∂µ i ∂ρ < 0 for any value of ρ.
We can use the Newton theorem that says that if :

1) f (ρ) is twice differentiable, V EXPERIMENTAL RESULTS
2) f (ρ * ) = 0, 3) ρ (0) is "close to" ρ * , 4) "f (ρ * ) = 0", then the series defined by [START_REF] Verdu | Minimax Robust Discrete-time Matched Filters[END_REF] converges to ρ * with a quadratic speed. It is easy to prove that I * = {1, ...p} (if the eigenvalues are sorted so that µ * 1 ≥ µ * 2 ≥ ... ≥ µ * N )). But we must be careful because that is not true at any step n for I (n) .

Then this algorithm converges to the solution ρ * of our problem.

C. Uniqueness of the solution

If we denote λ 1 the largest eigenvalue of B -1 A, then 1 ≤ ρ * ≤ λ 1 . In particular cases, it may be possible to find several subspaces of dimension p for which the SNR ρ is maximal; in fact, this is not a problem. In such a case, we can take an interest in finding a subspace of higher dimension than p with the same SNR, or we can add a new criterion to choose among those subspaces.

D. Practical remark

At each step, the identification of I (n) requires a heavy calculation. C N p different combinations have to be tested, which can quickly increase to an unacceptable number of calculations. The algorithm proposed can be improved significantly. Instead of searching the optimal set of eigenvectors in a systematic way, one can advantageously rearrange the eigenvalues of A -ρ (n) B in decreasing order:

µ (n) 1 > µ (n) 2 > ... > µ (n)
N (note that those values can as well be positive or negative), and take, at step n, I (n) = {1, 2, ..., p} (we have seen that I * = {1, 2, ..., p}): equation ( 14) becomes

ρ (n+1) = p i=1 t (n) i t At (n) i p i=1 t (n) i t Bt (n) i
. In theory, there is no reason for ρ (n+1) to be maximal, but in practice it so happens that ρ (n+1) reaches almost systematically the maximal value, and if not, reaches a value very close to it. In the neighbourhood of the solution, the convergence is assured. Such a change of the algorithm obviously highly decreases the sum of calculation. In terms of convergence to ρ * , there is a slight drop in the speed of convergence in terms of number of iterations. Globally, however this method converges to the solution and decreases the sum of calculation in a very significant proportion. To give a precise idea of the gain, for N = 21 and p = 5, C p N = 20349. The convergence of this modified algorithm has not been proved.

E. Conclusion

In this section, the convergence of the algorithm proposed has been proved.

For a given value of p, we have initialized the algorithm with ρ (0) = λ 1 the largest eigenvalue of B -1 A, saying it was reasonable to choose this value. What we are searching for is, p being fixed, the global maximum w.r.t. ρ (there exist other local maxima), and this maximum is necessarily the nearest from ρ = λ 1 (1 ≤ ρ * ≤ λ 1 ).

Obviously, an initialization of the algorithm with any value ρ (0) = λ 1 can lead the algorithm to find a local maximum.

V. Experimental results

Let us apply the CSMF method to the example described in Section III-F. Results are shown in Figure 3: the quality of the ROC curves increases from p = 1 to p = 4 (or p = 5 which gives more or less the same results than p = 4) and decreases from p = 6. The optimal result is obtained for p = 4 or 5. 1 , but also its superiority in comparison with the ESMF method. Now the optimal subspace E * p (here p = 4) has been found with the CSMF method which is a projection (a reduced-rank method). Nevertheless we did not use all the possibilities of classical filtering, and among all the basis of E * p , we can choose one with interesting properties after a linear filtering. For example, the ESMF gives preferential treatment to directions (spanned by the u i ) with the best SNR (λ i ) (like the Wiener filter): after this linear filtering, the power of noise is one in any direction u i .

Thus, we calculate the ROC curves obtained with the equation ( 2) where p is the subspace dimension so that all the basis vectors of E * 4 are taken into account. Results are shown in figure 5. A noteworthy improvement can be observed by using a simultaneous diagonalization technique in the optimal subspace calculated beforehand. We finally used a projection (to find the optimal subspace) and a linear filtering operation to improve again the detection.

VI. Conclusion

The method proposed in this paper takes its place in the set of methods of decomposition of signals on appropriate basis but also in methods.

When trying to detect stochastic signals with known covariance matrices but with no a priori knowledge on their probability density function, people usually try to project on the signal subspace (SVD,...). It is possible to take into account the structure (covariance) of the embedding noise: the SMF is used in such a point of view and in this case, a projection onto a p-dimensional subspace is made. In fact, this method is proven to be equivalent to a two-stage method: the whitening of the noise followed by the maximization of the SNR in a p-dimensional subspace. This method comes down to projecting the observation onto a subspace of dimension greater than one. However, there is no guarantee that the signal-to-noise ratio is maximum in the subspace spanned by these vectors.

In this paper, we calculate a subspace whose dimension is chosen a priori and which is optimal in the sense that the SNR ratio is maximized within. We prove, through theoretical examples, that this subspace is not necessarily those spanned by the vectors calculated by the SMF. Through ROC curves, practical experimentations illustrate the interest of such an approach.

We have shown, with a practical example, that a noteworthy improvement can be reached with the ESMF applied in the optimal subspace calculated beforehand by the CSMF. It confirms that such a method is an interesting and powerful one to perform detection.

Prospects of applications of the CSMF can easily be imagined in image processing or stochastic transient signals detection (like acoustic signals). An extension to the classification problem is possible. Of course, as this method is a reduced-rank one performing SNR maximization, it can be used for data compression or estimation and filtering. Thus, the CSMF was successfully used with real signals:

• analysis of sequences of IR images (SATIR) to qualify high heat flux components (carbon bricks used in the ITER project with the CEA Cadarache): detection of defects and classification of components [START_REF] Cismondi | Analysis of SATIR Test for the Qualification of High Heat Flux Components: Defect Detection and Classification by Signal-to-Noise Ratio Maximization[END_REF],

• detection and classification of textured images (like expanded polystyrene and textured paper or textured stones pictures for example): lot of images are texture images (forests, farming areas, ...) [START_REF] Borloz | Estimation, détection et classification par maximisation du rapport signal-à-bruit : le filtre adapté stochastique sous contrainte[END_REF],

• detection and localization of very high energy neutrinos by a passive underwater acoustic telescope (ANTARES European project) [START_REF] Juennard | Xerri Detection and Localization of Very High Energy Neutrinos by a Passive Underwater Acoustic Telescope Passive'08[END_REF],

• estimation of the sources in a specific blind source separation problem [START_REF] Xerri | Borloz Detection by SNR maximization: application to the blind source separation problem, 5th International Conference On Independant Component Analysis and Blind Signal Separation[END_REF].

Reduced-rank estimators and filters are important for a wide range of signal processing applications, among others when data or model reduction, robustness against noise or model errors is desired. This concerns known methods like the reduced-rank Wiener filter (RRWF) by Scharf [START_REF] Scharf | Matched and Adaptive Subspace Detectors When Interference Dominates Noise[END_REF], the reduced-rank maximum likelihood estimation (RRMLE) by Stoica-Viberg [START_REF] Stoica | Maximum Likelihood Parameter and Rank Estimation in Reduced-Rank Multivariate Linear Regressions[END_REF], the relative Karhunen-Loeve transform (RKLT) by Yamashita-Ogawa [START_REF] Yamashita | Relative Karhunen-Loeve transform[END_REF] or the generalized Karhunen-Loeve transform (GKLT) by Hua-Liu [START_REF] Hua | Generalized Karhunen-Loeve Transform[END_REF], which is used for data compression and filtering and is in fact nothing else than the RRWF also called low-rank Wiener filter in [START_REF] Scharf | Statistical Signal Proces.: Detection, Estimation and Time Series Analysis Reading MA[END_REF] section 8.4.

The choice of the optimal dimension p of the subspace of projection is an important question which must be examined in more detail in the future.

Fig. 1 .

 1 Fig. 1. ROC curves (input SNR=-14 dB) for different values of p

Fig. 2 .

 2 Fig. 2. µ i (ρ) for example 1 (N = 3)

Fig. 3 .Figure 4 Fig. 4 .

 344 Fig. 3. CSMF : ROC curves for different values of p

Fig. 5 .

 5 Fig. 5. ROC curves for SNR=-14 dB

  We denote U and ∆ λ the matrices such that AU = BU∆ λ : ). In this example, since t * 2 = u 2 , E * 2 is spanned by two eigenvectors of B -1 A, namely u 1 and u 3 which are clearly not those associated with the two largest λ i .Denote U and ∆ λ the matrices satisfying AU = BU∆ λ . = 1051, λ 2 = 0.45 and λ 3 = 0.01. The SNR obtained in E † 2 spanned by the eigenvectors associated with the two largest eigenvalues of B -1 A is ρ † = 149.64.

	00 0.80 0.64 0.80 1.00 0.80 0.64 0.80 1.00 1.00 0.60 0.36 0.60 1.00 0.60 0.36 0.60 1.00 0.6611 -0.7071 0.4170   , B = 1 3     . U =   0.3549 0.0000 -0.8076 0.6611 0.7071 0.4170 and λ 1 = 1.23, λ 2 = 0.56 and λ 3 = 0.46. The SNR   , obtained in E † 2 spanned by the eigenvectors associated with the two largest eigenvalues of B -1 A (u 1 and u 2 ) is ρ † = 1.06. The SNR obtained in E * 2 is ρ * = 1.118. If we denote (A -ρ   0.3370 -0.7071 0.6216 -0.8791 0.0000 0.4766 0.3370 0.7071 0.6216   , and µ Example 2 Let us consider two covariance matrices of non stationary processes. A = 1 3   0.0379 0.0379 0.1514 0.0379 0.0473 0.2650 0.1514 0.2650 2.9148   , and B = 1 3   1.2872 1.4658 0.1313 1.4658 1.6865 0.1629 0.1313 0.1629 0.0263   . U =   0.5219 0.7555 -0.3972 -0.5189 -0.6549 0.9156 0.6770 0.0192 -0.0622   , and λ 1   0.7111 -0.2531 -0.6559 -0.5845 0.3057 -0.7516 -0.3908 -0.9179 -0.0695   . and µ * 1 = -0.54, µ * 2 = 0.54 and µ * 3 = -150.05. E * 2 is spanned by t * 1 and t * 2 (we verify that µ * 1 + µ * 2 = 0 and µ * 3 = 1 -ρ 2 , a fortiori the eigenvector associated with the largest eigen-value of B -1 A which generates E * 1 . This example proves that a recursive algorithm w.r.t. p, that would calculate E * 1 and then E * 2 , etc... is not realistic. Conclusion From these simple examples, we immediately see that the optimal subspace E * p is not necessarily spanned by eigenvectors of B -1 A, and even when this is the case, the eigenvectors are not necessarily those associated with the largest eigenvalues. It is not possible to deduce E * p from E † p . What's more, it is not possible to find a recursive formulation on p to find E * p from E * p-1 : for example, the relationship E * 1 ⊂ E * 2 is not necessarily verified. Then we have to propose an algorithm to determine E *

* B) T * = T * ∆ * µ , then:

T * = * 1 = -0.0725, µ * 2 = -0.1186 and µ * 3 = 0.0725. E * 2

is spanned by t * 1 and t * 3 (we verify that µ * 1 + µ * 3 = 0 and

µ * 2 = 1 -ρ * The SNR obtained in E * 2 is ρ * = 151.05. We denote (A -ρ * B) T * = T * ∆ * µ , then: T * = * ). No eigenvector of B -1 A is contained in E *