Jihene Tounsi

Julien Boissière

Georges Habchi

Van-Dat Cung

A generic knowledge Model for SME Supply Chain Based on Multiagent Paradigm

.

According to some field investigations, three major features of the supply chain which integrates SME clusters arise. Firstly, a supply chain is a complex system, in particular in the SME context. This complexity is due to the number of autonomous actors and the number of SME networks which are linked and work together to achieve given processes and goals. Secondly, the SME are not often located in the same geographical area as it could be in a more classical supply chain. Finally, they face a lack of visibility over the entire supply chain as a result of the two previous characteristics. Indeed, sites only have local visibility but are coordinated with other sites through the flow of products. Due to this decentralized organization and limited view over the overall supply chain, studying the structure and the behaviour of the SME supply chain is a challenging task and even more so if sustainable considerations have to be taken into account.

Besides, the study of and the experimentation with the overall supply chain integration of SME clusters are difficult to implement on actual industrial systems without heavy investments from all the actors of the supply chain. Thus, in order to facilitate the analysis of the supply chain network, it is necessary to propose a

modelling solution which reflects the actual system and is able to simulate its behaviour. In light of this perspective, this chapter proposes a knowledge-model based on the Architecture Model Driven Engineering (ArchMDE [START_REF] Azaiez | Approche Dirigée par les modèles pour le développement de systèmes multi-agents[END_REF]) development process that aims to identify and model the domain concepts using the multiagent system. Hence, the work described here is a combination of two research areas. The first one (industrial engineering scope) proposes a modelling approach using different layers that represent different views of the system (i.e. the system refers to a supply chain). The representation of the domain concepts within the models allows one to capitalize on the know-how and then facilitates the reuse of the supply chain concepts in different contexts. The second one (the computer engineering scope) outlines the transition from the identification step of the domain structure to the study step for the dynamics behaviour of the domain concepts. This work aims to combine the domain concepts with the multiagent ones.

In this chapter, we highlight the research work through these two areas. Hence, the contents of this chapter are organized around four main issues:

 How can we build a reusable model?  Which methodology to adopt for domain conceptualization and what are the concepts?

 What are the steps to follow in order to agentify the domain metamodel?  How to move from the modelling stage to the implementation? To answer these questions, the first section introduces the ArchMDE development process and its contribution to this work. In the second section, we present the different steps applied to generate the conceptual metamodel. The third section highlights the agentification process and the dynamic behaviour integration of the different agentified concepts based on multiagent theory. Finally, the last section describes the transition from the modelling phase to the implementation one.

How to Build a Reusable Knowledge Model?

The heart of this research work is the modelling of SMEs supply chain by using multi-agent systems in order to build a reusable, flexible and secure knowledge model. To reach this aim, advances in the field of computer engineering, especially those dealing with multi-agent paradigm appear to be a promising approach.

To this end, we have adopted a modelling approach called ArchMDE and proposed a PhD thesis [START_REF] Azaiez | Approche Dirigée par les modèles pour le développement de systèmes multi-agents[END_REF]. This approach is based on Model-Driven Engineering (MDE) [START_REF] Kent | Model-driven Engineering[END_REF] which founds its developing process on producing several interrelated models. In the ArchMDE approach, two types of metamodels are identified: a domain metamodel that describes functional concepts and properties related to a particular domain (e.g. a SME supply chain) and a computer modelling metamodel (e.g. a multiagent system). The combination of both metamodels will generate an agentified metamodel, that constitutes the starting point of the conceptual models. From this last metamodel, different functional models are de-scribed in order to introduce the functionalities of the system (Fig 1). Finally, the use of a platform metamodel is necessary to generate the program code. [START_REF] Tounsi | A multiagent approach for modelling SMEs mechatronic supply chains[END_REF] This approach is of great interest to fill in the existing gap between the design and the implementation phases. The following sections describe the step-by-step approach through those main axes.

 The methodology to generate the domain metamodel.



The agentification process which is introduced through the multiagent metamodel and the analogy between the multiagent concepts and the domain ones. The agentification is achieved by integrating the dynamic behaviour into the agentified SME supply chain metamodel.

 The implementation phase focuses on the transition from the modelling step to the encoding one.

Conceptual Domain Metamodel

According to the ArchMDE development process, the first modelling step involves the definition of the domain conceptual model. This step leads to the identification of the main concepts of SME supply chains. To achieve this objective, we follow a methodology based on existing conceptual modelling visions in the literature [START_REF] Tounsi | A conceptual model for SME Mechatronics supply chain[END_REF]. In this methodology, the visions are organized into three steps. Each step addresses concepts related to supply chains. These concepts and their relationships will then be gathered within a domain metamodel that will be expressed using Unified Modelling Language (UML). The following section presents this methodology.

Conceptual modelling methodology

To identify the properties and concepts of the supply chain domain, an incremental methodology combining three visions is proposed: product vision, structure vision and process vision. In each step, a vision is applied to build or to refine the conceptual model. The result of each step (intermediate model) is the input of the next one. Therefore, at the end of the three steps, a final architecture of the conceptual model is generated (Fig 2).

Fig. 2

Conceptual modelling methodology framework [START_REF] Tounsi | A conceptual model for SME Mechatronics supply chain[END_REF] Step 1: Product Vision This vision considers the supply chain dedicated to a particular product (or a family of products) from the raw materials through to the final goods. It focuses on the product flow to define the environment and organizations involved in its management [START_REF] Thierry | Gestion des chaînes logistiques : Modèle et mise en oeuvre pour l'aide à la décision à moyen terme[END_REF]. In the methodology framework, the Product Vision leads to the construction of a first abstract model of the supply chain involving the environment and organizations:  The environment is characterized by the physical flows and the different steps of the product transformation as well as the related disturbances.



The organizations are the entities carrying out one or several product transformation stages and the physical flow management. The organizations involved can be a network of firms that collaborate to accomplish one or several transformation stages.

Step 2: Structure Vision This vision has been proposed by [START_REF] Cooper | Supply chain management: more than a new name for logistics[END_REF]. It considers the architecture of the supply chain, made up of: actors (decision-making actors and synchronization actors), network structure (roles in the network and the number of actors for each role) and relationship characteristics between actors.

Structure Vision

-Defines the environment architecture -Defines the organization's architecture

Intermediary conceptual model

Product Vision

-Defines the environment type -Defines the organisation types

Abstract model

Process Vision

-Identifies the processes -Integrates the processes into the model

Conceptual Model

Thus, on the basis of the abstract model provided by the previous step, the Structure Vision details the organizations involved and the physical environment:  The environment is the part containing the physical flow. Therefore, the product flow and the resources used to achieve its transformation have to be described.

 The organization consists in identifying and prioritizing the actors in the network according to their involvement in the different levels of decision-making as well as the tasks they will be assigned. The information flow management that depends on the decision-making level is also considered in this approach.

At this step, a more detailed intermediate model is built.

Step 3: Process Vision This vision is based on the process classification according to the decisionmaking levels [START_REF] Stevens | Integrating the supply chain[END_REF][START_REF] Chopra | Supply Chain Management: strategy planning and operation[END_REF]: strategic, tactical and operational.

While applying Process Vision, the various categories of processes are identified and integrated into the previous intermediate model. This can be done according to the decision level but it also depends on the relationships between the actors. These relationships can be classified into two categories:  Synchronization: contains processes for exchanging information and physical flows according to a process program developed and predefined by the decision-making layers.

 Management and control: contain processes that ensure suitable decision implementation in the perspective of a continuous improvement of processes in terms of added value.

This step leads to a refined conceptual model of the supply chain.

Domain Model Concepts

This section presents the concepts that constitute the domain model. By applying the methodology described in the preceding section, several concepts, processes and the architecture of the model were identified. Based on these concepts, a metamodel of supply chain is proposed.

Step 1: Applying Product Vision By applying the Product Vision, a first abstract model of the supply chain is built. It is composed of (Fig 3 ):  Environment: represent the part allocated to the product flow and management through the internal resources as well as the external elements able to influence supply chain activities.

 Sub Supply Chain (SSC): represents a group of SME which collaborate to achieve an internal aim and/or the overall objective of the supply chain. The SSC is responsible for the management of the product flow in a certain stage of its life cycle.



Perimeter of influence: represents the visible part of the environment to the SSC on which it can act by internal conferring (if the action does not disturb the environment located outside its visibility) or by conferring with another SSC.

 Shared perimeter of influence: represents the area of the flow transfer between two SSCs. It is a shared zone where SSC coordinates their activities to allow the flow transfer. [START_REF] Tounsi | A multiagent approach for modelling SMEs mechatronic supply chains[END_REF] Step 2: Applying Structure Vision By applying Structure Vision, the previous abstract model is refined. The internal architecture of the SSC and the visible part of the environment (the perimeter of influence) are described. As showed in Fig 5, the SSC model and its environment are based on three layers representing the different decision-making levels.

SSC( SSC(

Interac-

Perimeter of influence

Perimeter of influence

Shared perimeter of influence

Environment

Visibility Visibility

Fig. 5 Layers of the SSC [START_REF] Tounsi | A conceptual model for SME Mechatronics supply chain[END_REF] Each layer involves particular concepts and plays a specific role in the SSC:  The Monitoring System is the intelligent layer of the SSC. It controls and monitors the two other layers through the information provided by the Execution System. Monitoring Actors (MA) modelling the intelligent actors of SSC are the main elements of this layer. They establish metrics to evaluate the performance of the group and consequently act on the other two layers. Hence, MA are the components responsible of controlling and decision-making into a SSC and of the coordination of the activities for the overall supply chain.

 The Execution System is the reactive layer of the SSC. It has two main roles: (i) it ensures the synchronization of the physical flow according to the information gathered from the Physical System, (ii) it observes and corrects the Physical System if a disturbance occurs. In abnormal situations, the Execution System refers to the Monitoring System for coordination and decision-making. Executive Actors (EA) are the main entities of this layer. An EA mainly models the reactive actor in the Execution System. However, occasionally MA can appear in this layer with reactive behaviour.

 The Physical System is the visible part of the SSC environment. It corresponds to the SSC's perimeter of influence. This layer is composed of nondecisional elements controlled by the other two layers of the SSC. Two main concepts are identified: the Moving Entity (ME) modelling the product flow and the Resource modelling production means.

Fig 6 shows the first conceptual abstract model refined in a domain metamodel. On one hand, we have integrated the identified concepts of each layer. On the other hand, an abstract class "Actor" is added to the metamodel for implementation purposes. Indeed, the "Actor" class defines the structural characteristics and behaviour of a decisional entity. Thus, both the EA and the MA inherit from this class. However, the "EA" class defines the specific characteristics of an executive actor and likewise for the "MA" class.

Measurement

..n 1 1..n 1 SC SSC 1..n 1 1..n 1 1..n 1 1..n 1 1 1 1 1 1 1 1 1 Environment 1 1 1 1 Resource ME Physical System 1 1 1 1 1..n 1 1..n 1 1..n 1 1..n 1 1..n 1 1..n 1
Fig. 6 Intermediary domain metamodel [START_REF] Tounsi | A multiagent approach for modelling SMEs mechatronic supply chains[END_REF] Step 3: Applying Process Vision The object of the last step is to identify and integrate the different kinds of processes into the model. Table 1 gives a classification of the processes identified according to their role in decision-making. In the Physical System, the Physical Processes (PhP) have been identified. A PhP describes the sequence of the processing stages of a product. It is a concept to be integrated within a domain metamodel in order to define the tasks that can be handled by the Execution System.

-Define the transformation routings of products

The processes identified in both Monitoring and Execution Systems are management processes. Hence, they represent the dynamic behaviour of the SSC. This behaviour is induced by control and monitoring decisions that come from either the SSC or the overall supply chain. This behaviour is basically a communication mechanism (coordination, collaboration or cooperation).

In order to model management processes and communication mechanisms, more informational elements are needed for EA and MA to ensure their role in the domain model. Thus, decisional actors of the SSC (EA and MA) need three conceptual elements that consolidate their internal architecture:

 Indicator: is used by actors for two different tasks. Indeed, the EA control and detect Physical System deviation by comparing the value of an indicator with its fixed objective. As for MA, they evaluate the internal performance of the SSC but also in the overall supply chain.

 Action: Actors apply actions when facing indicator deviation.



Organizational Knowledge: is an actor's database that stores information about his acquaintances. For example, if an actor "A" is an acquaintance of an actor "B" this means that "B" owns information about the identity, the behaviour, the capabilities and the resources of the actor "A". Reciprocally, the actor "A" owns the same information about the actor "B". According to this, each actor (EA and MA) owns knowledge about resources of all actors in the same SSC. However, the MA involved in the overall supply chain have additional internal acquaintances, each MA owns limited knowledge about the other MA of the overall supply chain. Note that this knowledge requires continuous updating.

In the same way, the intelligent behaviour of the MA requires the definition of other conceptual components:

 Objective: models the strategic goal of the SSC. According to this aim, the SSC coordinates its activities with other SSCs in the overall supply chain.

 Constraint: is a knowledge that an actor must consider to reach the goal of the overall supply chain or the SSC's one.

Through the Process Vision, the previous metamodel and its concepts are refined by integrating identified concepts. Figure 7 [START_REF] Tounsi | A multiagent approach for modelling SMEs mechatronic supply chains[END_REF]

1 1 coordinate SSC 1..n 1 1..n 1 1 1 1 1 1..n 1 1..n 1 1..n 1 1..n 1 1 1 1 1 1 1 1 1 1..n 1 1..n 1 monitors 1..n 1..n 1..n 1..n 1..n 1..n 1..n 1..n Fig 7. Final domain metamodel

Agentified SME supply chain metamodel

In this section, the domain metamodel is merged with an agent metamodel using the ArchMDE methodology [START_REF] Azaiez | Approche Dirigée par les modèles pour le développement de systèmes multi-agents[END_REF]. Thus, this section is divided into two parts: the first one outlines the properties of each multiagent concept according to the Agent-Interaction-Environment-Organization (AIEO) approach and the second part highlights the agentification of the domain metamodel.

Multiagent metamodel

The AIEO approach breaks the whole multiagent system down into four views: Agent view, Interaction view Environment view and Organization view. Fig. 8 Agent metamodel [START_REF] Azaiez | Approche Dirigée par les modèles pour le développement de systèmes multi-agents[END_REF] The Agent view defines the agent metamodel composed of the following concepts:

 "Agent" identifies different kinds of agent according to the decisionmaking capacity of the agent (reactive agent, cognitive agent and hybrid agent).

 "Cognitive agent" defines an agent with cognitive abilities. The metamodel highlights the main concepts modelling the BDI agent (Belief, Desire, In-tention=plan).

 "Reactive agent" defines an agent with reactive abilities to respond to unpredictable events.

 "Hybrid agent" defines an agent with hybrid intelligence (cognitive and reactive abilities).



"Goal" defines the aim that an agent should achieve.  "Knowledge" and "Norm" define all the knowledge and norm necessary for the agent to achieve its goal.

 "Plan" represents an action plan implemented by the agent. The plan is composed of one or several elementary actions.

 "Reactive action" is an action implemented by the reactive agent. The Interaction view describes the dynamic relations between the agents. This interaction is a structured exchange of messages between the agents through a specific protocol or language. Thus the interaction metamodel highlights the following concepts:

 "Interaction protocol" represents the interaction protocol adopted by the agents.

 "Communicative action" represents an elementary action of communication that is part of the "interaction protocol".

 "Message" is a set of information exchanged between the agents through the "interaction protocol". The agent interprets the message based on the communicative action.

The Environment view focuses on all the elements external to the agent allowing it to reach its goal or activate its behaviour through events. The elements belonging to the environment metamodel are as follows:

 "Active resource" represents the resources that activate the behaviour of the agent by generating events or triggers.

 "Passive resource" defines the resources the agent needs to accomplish its task.

 "Event triggered action" represents events that resources activate. An event is composed of one or more tasks.

Finally, the Organization view describes the structure of the whole system. The organization metamodel is made up of the following concepts:

 "Organization" defines the system topology (hierarchy, group or market).  "Role" represents different roles that the agent could play.

Agentification of the domain metamodel

This step of the ArchMDE methodology consists in merging the multiagent metamodel with the domain metamodel. Hence, on the one hand, a metamodel defines a multiagent system according to the Vowel approach (Fig 8). On the other hand, a domain metamodel describes the supply chain in SME context (Fig 7). A correspondence between the multiagent concepts and those of the domain is then carried out according to their properties and their roles in the metamodel. Table 2 summarizes the correspondence between these concepts in order to achieve the agentified metamodel for SME supply chains. According to the collected information and the history of the situation and action, the group of MA monitors the SSC to reach a goal and accomplish its activity. Objective Goal Desire

A SSC has a goal to reach. This goal is coordinated with other nodes' goals. In addition, each MA has a personal aim for each indicator. This kind of "Objective" is modelled by the "Desire" of the BDI agent (MA).

Indicator Belief Perception

The agents act on the environment according to the indicator measures. In this case, an indicator is modelled by the "Perception" of an agent.

However, a MA monitors the SSC according to the history of these measures. So, an "Indicator" is modelled by the "Beliefs" of BDI agent. Action Plan It is an action or a set of actions to apply when facing a disturbance. Knowledge Knowledge It is all the knowledge needed by the agents to act in an appropriate way.

Organizational knowledge

Knowledge

Each agent has a list containing the information about other agents from the same SSC or the overall SC. This list stores knowledge about the name of the agent, the task that it performs and its resources.

Constraint Knowledge

The MA make decisions according to their objectives and their beliefs. At the same time, there are some constraints (about product or other SCs where the group is involved) that the group of MA must take into account when making decisions.

After the agentification process, we obtain an agentified domain metamodel as presented in Fig 9 . The domain metamodel and the multiagent one are separated for more clarity.

The integration of processes into the metamodel

Up to now, the static part of the domain metamodel has been created. In this section, we define the dynamic behaviour of the concepts based on the multiagent tools and theory. Indeed, this dynamic is described by the implementation of interaction protocols according to the process vision and the communication mechanisms.

Firstly, the process vision allows us to define two scenarios: (1) the synchronization of the physical processes and (2) the monitoring and the control of processes.

Secondly, a communication mechanism is "a framework formalizing interaction between different actors in the network according to their managerial relationship characteristics" [START_REF] Tounsi | A conceptual model for SMEs mechatronics supply chain[END_REF]. The study of the domain identifies two kinds of communication framework. Indeed, in the overall supply chain, SSC coordinate their activities in order to achieve the common objective of the overall supply chain. Within the SSC, the actors collaborate to achieve a local goal. This section describes the different protocols implemented in the agentified domain metamodel taking into account the different scenarios of the process vision and the communication mechanisms.

Fig. 9

The agentified SME supply chain metamodel [START_REF] Tounsi | A multiagent approach for modelling SMEs mechatronic supply chains[END_REF]

Synchronization of the physical processes

The SSC is responsible for the synchronization of the Physical System involved to achieve its task. This activity consists in applying a communication protocol relative to the nature of the interaction framework. In this section, the collaboration and coordination processes are described in order to be implemented in the Execution System and Monitoring System, and then to synchronize the Physical System.

Integration into the execution system

In accordance with the agentified domain metamodel, the Execution System is responsible for the synchronization of the physical process (PhP) in common situations. Indeed, Executive Actors (EA) which are reactive agents, synchronize PhP by taking into account the availability of resources. The following sequence highlights the EA's behaviour:  If the request is a negative response for the collaborative demand that the EA initiated, the EA sends a request to the coordinator agent (MA coordinator) of the monitoring system.

 If the request is a synchronization request coming from the Moving Entity (ME) or a collaborative request coming from another initiator agent, the EA checks the availability of the resources concerned.

 If the resource is available, the EA carries out its task, updates the state of the ME and informs the other agents within the executive system and the coordinator agent at the end of the action.

 If the resource is unavailable and the EA has been solicited by another executive agent to achieve the task, it sends a failure request to the initiator.



If the resource is unavailable and the EA is in charge of the task, then it seeks within its organizational knowledge an agent from the Execution System of the SSC that might have the necessary resource.

 If the agent finds another agent within its organizational knowledge that can handle the task, it delegates the responsibility of the task. In this case, the collaboration process of the agent concerned will be activated and follows the same sequence.

 If the agent does not find another agent who has the resource necessary to handle the task, it sends a request to the coordinator agent. This agent is a monitoring agent (MA) that receives requests from the Execution System. The MA sends the information to other monitoring agents in the SSC in order to find a solution.

Integration into the monitoring system

In unusual situations1 , the Executive System refers to the Monitoring System. In this case, the group of MA evaluates the situation according to the defined objective and establishes an action plan. If the objective is not reached, the MA needs to consult other SSCs to find a suitable solution. Thus, the agents adopt the protocol based on Contract Net Protocol to provide the coordination of the objectives. The synchronization protocol can be described according to the following steps:  In the Monitoring System, a monitoring actor is responsible for checking all the requests received and sending them to other MA in the layer. Three kinds of requests can be distinguished: (1) EA request, (2) reply to a help request or (3) a help request from another SSC in the overall supply chain.

 If a MA coordinator receives an EA request then it sends the information to other MA. In this case, the group evaluates the situation according to the SSC's objectives. Two cases may arise: (1) the problem has no impact on satisfying the SSC's objectives or (2) the objective is deviated.

 If there is no impact on the objective, the MA tries to find an internal solution according to its desire, belief and constraints. If a solution can be found, the MA coordinator sends the actions plan to the Executive System.

 If the objective is deviated or an internal solution cannot be reached, the group of MA sends a Help Request to other SSCs via the MA coordinator and waits for the responses.

 When all the replies have been received, the MA coordinator ranks them according to the reception date of the request. The list of responses will be sent throughout the Monitoring System. According to their beliefs, desires and constraints, the group of MA chooses the most suitable answer and diffuses the action plan to the Execution System. In this case, the EA updates the state of the ME.



If the request is a help request from another SSC, the MA coordinator sends the request to the Monitoring System. In this case, the Monitoring System evaluates the demands according to internal criteria (objective, constraint, belief, and desire). If the SSC can provide assistance, it makes an offer to the SSC initiator or it sends a negative response.



The SSC initiator chooses the suitable offer and sends a confirmation to the selected SSC and a cancellation response to other bids.

The following figure (Fig 11) shows the sequence of messages between the SSCs. This diagram represents the coordination process in the overall supply chain in order to synchronize the physical flow in the case of a disruptive case (SSC cannot reach the internal aim).

Monitoring and control protocol

This protocol describes the conditional preventive (based on measurements) and corrective (in case of disturbance) monitoring and control in the SSC or in the overall supply chain. The monitoring and control protocol is based on perform-ance evaluation in both layers of the SSC: the Execution System and the Monitoring System.

At the end of the synchronization protocol, each actor updates the indicator measurement by evaluating its activity and related resources. In addition, the Monitoring System evaluates the local activity of the SSC and participates in the improvement of the performance of the overall supply chain. Thus, the following sequence describes the monitoring and control protocol in the SSC's layers:

 At the end of its synchronization task, the EA evaluates the performance of its activity and related resources (the allocated space of the environment to the EA).

 According to this perception, the EA refers to the indicator base in order to detect a disturbance.

 If the EA finds a deviation, it seeks the cause of the disturbance.



If the deviation is a common situation, the EA selects the appropriate action plan to solve the problem and applies it to the environment. After that, it sends measurement (or perception) to the MA coordinator which, in turn, sends the information throughout the Monitoring System. Then, each MA updates its belief.



If a new situation occurs, the EA sends a failure control message to the MA coordinator. Then, the MA sends the information throughout the Monitoring System and each MA updates its belief.

 In this case, the Monitoring System analyzes the situation according to internal criteria (beliefs, desires, objectives and constraints).

 If the problem needs corrective maintenance, the Monitoring System generates an action plan and forwards it to the Execution System. The actors of the SSC update their bases of actions.

 If the disruption does not affect the SSC, the Monitoring System applies a preventive action plan to avoid future disturbance.

Fig Fig. 12 Sequence diagram of the monitoring and control protocol [START_REF] Tounsi | A multi-agent knowledge model for SMEs mechatronic supply chains[END_REF]

The Implementation Phase

The last phase of the ArchMDE development process deals with the transition from the modelling phase to the implementation one. This last phase is mainly divided into two steps. The first one concerns the refinement of the agentified domain metamodel by integrating the information necessary for the implementation. The second step introduces the choice of the development platform. In the next section, we focus on the two steps of the implementation phase of the ArchMDE development process.

Refinement of the Agentified Domain Metamodel

The steps of integration protocols identify the abstract patterns that describe the dynamic behaviour of the agentified domain concepts. The refinement consists of the integration of the attributes and the methods that define the architectural and the behavioural proprieties of each concept. This step provides a final class diagram, the Implementation Metamodel. The encoding of this last result leads to a dedicated simulation platform.

In the implementation metamodel, we have refined the domain metamodel and integrated multiagent concepts that are essential to the agent's behaviour. The EA (reactive agent) and MA (cognitive agent) inherit from Actor class the common characteristics. Nevertheless, each one implements the method run() that describes his behaviour and the method HandleMessage() that allow him to read and to construe the received message.

The choice of the development platform

This section focuses on the choice of a development platform according to several objectives of the implementation. Indeed, obtaining the final pattern (the implementation metamodel) completes the modelling process. At this level, the behavioural and architectural characteristics of the proposed concepts are implemented using a development platform. The most important modelling criteria that the implementation platform must allow are as follows: (i) the implementation of cognitive agents (MA) and reactive agents (EA), (ii) the communication between agents (sending messages), (iii) the integration of the agents group (Monitoring System and Execution System) and the modelling of the external environment (Physical System), and (iv) a graphical user interface for easy setup and the observation of the simulation results.

According to these criteria, we have studied and compared two multiagent platforms: Jade (http://jade.cselt.it) and Madkit (http://www.madkit.net).

Jade (Java Agent Development Platform) is a multiagent platform fully encoded in JAVA. It allows the modelling of agents based on predefined patterns communicating through messages. The Jade software simplifies the implementation through a graphical user interface (remote GUI). However, it does not allow the implementation of a group of agents which is a major modelling criterion for our metamodel. Consequently, Jade does not correspond to our specifications.

Madkit is a modular and scalable multiagent platform also written in JAVA. The main reasons for taking an interest in this software are that it: (i) provides an API (Application Programming Interface) to enable the development of agents that communicate through sending messages, (ii) allows one to develop agents located in groups and play roles in the organization and (iii) offers a full set of facilities for launching, displaying, developing and monitoring agents and organizations [START_REF] Gutknecht | MadKit: Une expérience d'architecture de plate-forme multi-agents générique[END_REF]. However, Madkit does not allow one to draw the external environment as a set of objects (object in oriented-object programming theory). Indeed, each concept must be an agent in order to communicate into the application. Nevertheless, we can use the Madkit platform with a JAVA environment to integrate external classes. Thus, Madkit merged with JAVA environment can be considered as an implementation way. Due to lack of time and knowledge in this area, we have not investigated further, but this solution can be considered as a future perspective to implement the knowledge model.

Given our computer skills, we have decided to implement our own simulation platform using JAVA. In this platform, we have developed agents, groups of agents and structured a peer-to-peer communication between them. The elements external to the agent (i.e.: resources, indicator, action, PhP, etc.) are encoded as objects. The product (ME) is an active entity. It is encoded as an object that triggers events.

Conclusion

For reasons of complexity in terms of size, communication protocols and decision-making strategies with decentralized behaviours of actors, modelling and simulation of supply chains is a tricky task for researchers and supply chain managers. In this chapter, our objective is to propose a natural way for supply chain modelling while considering the organisational issues and the managerial relationships. The resulting knowledge model requires three main features: reusability, simplicity and genericity. Our main theoretical contributions are on this conceptual model and its building methodology.

To reach this goal, we have adopted ArchMDE as a development process. The first step of the development process focused on the identification of the major concepts. These concepts are consistent with both domains: the inherently decentralized supply chain and multi-agent systems. The high level of conceptualization of these concepts tends to provide genericity and reusability to the proposed model. In the second step, we have presented a mapping within the proposed concepts. This mapping translates the supply chain concepts into MAS concepts. This aims at simulating the dynamic behaviour of the supply chain based on multiagent tools. Major communication protocols between actors (or agents) are also drawn to emphasize the importance of different kinds of exchange mechanisms. To keep the model as generic as possible only a few generic protocols are presented here.

In terms of implementation, the concepts and the protocols presented in this chapter have been developed as a generic Supply Chain Multi-Agent System using JAVA. A prototype of a 2-echelon Supply Chain with one manufacturer and one logistician including (1) respectively production and inventory planning processes, and (2) two-actor negotiation and coordination processes has been tested with success [START_REF] Ogier | Supply Chain Performance in the case of Decentralized Planning[END_REF].

Our aim is to be able to simulate a real multi-actor Supply Chain with dynamic behaviour. Thus, more internal optimization processes and coordination protocols are under study.

Fig. 1

 1 Fig.1ArchMDE development process[START_REF] Tounsi | A multiagent approach for modelling SMEs mechatronic supply chains[END_REF]

Fig. 3 Fig. 4

 34 Fig.3The abstract model[START_REF] Tounsi | A conceptual model for SME Mechatronics supply chain[END_REF] Fig4shows the domain metamodel which reflects this conceptual abstract model using UML.

  presents a UML representation of the final domain metamodel for the supply chain in the SME context. It corresponds to the final conceptual model with its associated concepts regardless of computer technologies.

  Fig 8 shows multiagent concepts according to each view and the links between them.

  Fig 10 describes the Executive Actor's (EA) behaviour in its role of synchronizing the physical flow (ME). An EA receives a request and reacts according to its type. Three types of request can be distinguished: (1) a ME request, (2) a col-another EA or (3) a negative response to a collaborative action initiated by the agent itself.

Fig. 10

 10 Fig.10EA synchronization behaviour(Tounsi et al. 2009a)

Fig. 11

 11 Fig.11SSC coordination process[START_REF] Tounsi | A multiagent system for production synchronization in SME mechatronic supply chain[END_REF]

  Fig 13 represents the Implementation Metamodel. The class diagram focuses on the main methods and attributes. Each private attribute has accessor (get) and mutator (set) methods. However, these methods are not visible in order to lighten Fig 13. The Actor class defines the global architecture and behaviour of an agent.

Fig. 13

 13 Fig.13The Implementation Metamodel[START_REF] Tounsi | A multi-agent knowledge model for SMEs mechatronic supply chains[END_REF]

Table 1

 1

		Process Classification
	SSC Layer	Process Family	Role
	Monitor-	Strategic Processes (SP)	-Coordinate long term decisions
	ing Sys-	Monitoring and	-Monitor SSC activities
	tem	Control Processes	-Drive and evaluate SSC performance in
		(MCP)	the overall supply chain
	Execution	Operational Control	-Synchronize and control the physical sys-
	System	Processes (OCP)	tem
	Physical	Physical Processes
	System	(PhP)

Table 2

 2 Correspondence between domain and agent concepts

	Domain con-	Multiagent	Description
	cepts	concepts
	Supply Chain	Multi-Agent	By analogy, the root of the domain metamodel
	(SC)	System	corresponds to the root of the multiagent sys-
		(MAS)	tem.
	Environment	Environment	In both metamodels, the environment is the
			physical space defining all things that are exter-
			nal to the agents and necessary in order to man-
			age the SC.
	Sub Supply	Organization	It is an organization made up of two groups of
	Chain (SSC)		agents.
	Physical Sys-	Resource	It is all the resources needed for one agent or a
	tem		group of agents to manage the group (perimeter
			of influence).
	Resource	Passive Re-	It is a resource allocated to the agent to perform
		source	its task.
	Moving En-	Active Re-	The ME represents the product in circulation. It
	tity (ME)	source	activates the behaviour of the reactive agents.
	Physical	Task	It is a task or a physical activity to be handled
	Process (PhP)		by reactive agents.
	Monitoring	Group	It is a group of cognitive agents which collabo-
	System		rate in the SSC and coordinate the activity of
			the organization with other organizations.
	Execution	Group	It is a group of reactive agents which collabo-
	System		rate in the SSC.
	Actor	Agent	An actor can be a cognitive agent or a reactive
			agent according to its decisional characteristics.
	Executive	Reactive	EA perceives the physical system and acts on it
	Actor (EA)	Agent	according to the observation.
	Monitoring	Cognitive
	Actor (MA)	Agent

  12 shows the protocol for monitoring and control through an UML sequence diagram:

	: EA	: Indicator	: Action	: Environment MA _coordinator	: Monitoring System
				: MA
	t ask_end()
		perceives( )
		return (list_perception)
	consult(id_perception)
	return (alertmin, alertmax)
	compare(perception, alertmin, alertmax)
	select(id_indicator, pert urbation_cause)
	[action exist] retrun (action)
		[act ion ex ist ]ac t (act ion)
		s ends _Lis t(indicat ors , values)
				sends _List(indicators, values)
					Updates _Belief(id_indicator, value)
	[action doesn't exist ]send(failure_message,id_indicators, values)
				sends _List(indicators, values)
					Updates_Belief(id_indicator, value)
				analyze(Belief, Desir, Intention, Constraint)
				sends(action_plan, id_indicator, perturbation_cause)
		sends(action_plan, id_indicator, perturbation_cause)

Unusual situation is occurred when Executive System cannot propose a solution for a happened problem in the physical system.


Matthieu Landès
email: landes@ipgp.jussieu.fr

Fabien Hubans
email: fabien.hubans@obs.ujf-grenoble.fr

Nikolai M Shapiro
email: nshapiro@ipgp.jussieu.fr

Anne Paul
email: anne.paul@obs.ujf-grenoble.fr

Michel Campillo
email: michel.campillo@ujf-grenoble.fr

Origin of deep ocean microseims by using teleseismic 1 body waves 2

Recent studies of oceanic microseisms have concentrate on fundamentalmode surface waves. Extraction of fundamental-mode Rayleigh and Love wave Green functions from station-station correlations of ambient seismic noise has recently been demonstrated to be a very powerful tool for imaging of the Earth's crust and uppermost mantle.

In this study we concentrate on energetic arrivals in two frequency bands around the primary (14s) and the secondary (7s) microseismic peaks that appear at near-zero times in noise cross-correlations. Thanks to a polarisation analysis of data from the the ETSE network (Turkey), we identify this "near-zero time" signal as an upcoming P wave in the secondary microseismic frequency band (5-10s). In a second step, analysing noise cross-correlations from three different arrays ( in Yellowstone, in Turkey and in Kyrgyzstan), we determine the origin of these signals by means of beamforming analysis and its projection on the Earth.

Our results show that, in the 0.1-0.3 Hz frequency band, the energetic "nearzero" time arrivals in seismic noise cross-correlations are mainly formed by teleseismic P, PP, and PKP waves. Generation of this ambient body waves in the secondary microseismic band presents a marked seasonal behaviour with sources located in southern and northern oceans during summer and winter, respectively. Moreover, body wave array analysis is accurate enough to confirm that significant amount of the microseism energy is generated far from the coast in deep oceans.

Introduction

Recent years witnessed a strong interest in studying background seismic noise. One of the reasons for this interest is the possibility to extract deterministic Green functions from correlations of a random wavefield that can be proved mathematically with different approaches (e.g., [START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF]; [START_REF] Snieder | Extracting the Green's function from the correlation of coda waves: A derivation based on stationary phase[END_REF]; [START_REF] Gouédard | Crosscorrelation of random fields: mathematical approach and applications[END_REF]) and that has been demonstrated in acoustic laboratory experiments (e.g., [START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF]; Derode et al. [2003]). Application of this principle to large amount of continuous digital seismic records provided by modern networks provides us with new approaches for seismic tomography (e.g., [START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]; [START_REF] Shapiro | High-Resolution Surface-Wave Tomography from Ambient Seismic Noise[END_REF]; [START_REF] Sabra | Surface wave tomography from seismic ambient noise in Southern California[END_REF]; Yang et al. [2007]; [START_REF] Stehly | Cross-correlation of random fields: mathematical approach and applications[END_REF]) and monitoring (e.g., Sens-Schönfelder and [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia[END_REF]; Brenguier et al. [2008b]; Brenguier et al. [2008a]). Detailed analysis of high-quality continuous records also allows us to better understand the origin of the ambient seismic noise and its relation to oceanic and atmospheric processes. New methods of noise-based seismic imaging and monitoring are based on a principle that the Green function between two points can be extracted by correlating a random wavefield recorded by receivers located in these points. In other words, one of the two receivers can be considered as a virtual source recorded by the second receiver. This principle is especially attractive when applied in context of random wavefield recorded by a network of numerous recorders.

In this case, by computing all possible inter-station cross-correlations it becomes possible to place virtual sources at every receiver location and to have their records by the whole network resulting in a very dense path coverage. The deterministic waveforms (Green functions) extracted from the cross-correlations can be used then to measure and to in- The noise based Green function reconstruction implies, however, some strong hypothesis about the noise modal composition. A perfect reconstruction can be achieved for an ideally random and equipartitioned wavefield [START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF]; [START_REF] Sánchez-Sesma | Retrieval of the Green's Function from Cross Correlation: The Canonical Elastic Problem[END_REF]). In a case of seismic noise it would imply that its sources should be distributed randomly and homogeneously in volume. This is obviously not the case for the real ambient seismic noise within the Earth. First, most of its sources are located on the surface resulting in stronger presence of fundamental mode surface waves and their relatively easy reconstruction from inter-station cross-correlations. For the same reasons, extracting body wave Green functions from noise cross-correlations remains challenging.

Second, the distribution of noise sources is not perfectly random and homogeneous. Background seismic oscillations are mostly generated by the forcing from oceanic gravity and infragravity waves. The interaction between these oceanic waves and the solid Earth is governed by a complex non-linear mechanism (e.g., [START_REF] D R A F T ; D R A F T Longuet-Higgins | A Theory of the Origin of Microseisms[END_REF]) and, as a result, the noise excitation depends on many factors such as the intensity of the oceanic waves but also the intensity of their interferences as well as the seafloor topography (e.g., [START_REF] Kedar | The origin of deep ocean microseisms in the North Atlantic Ocean[END_REF]). Overall, the generation of seismic noise is strongly modulated by strong oceanic storms and therefore, has a clear seasonal and non-random pattern.

Distribution of noise sources homogenizes when considered over long times (more than one year). The homogenization and randomization of the noise wavefield is also enhanced by the scattering of the seismic waves on the small-scale heterogeneity within the Earth.

Also, because of the stationary phase principle, a cross-correlation of the noise recorded

D R A F T December 1, 2009, 5:06pm D R A F T X -5
by two receivers is dominated by contribution from sources located in vicinity of the line connecting these receivers. Therefore, even without a perfectly homogeneous distribution, a presence of sufficient amount of favorably located noise sources results in relatively high quality reconstruction of fundamental mode surface waves. As a consequence, reconstructing surface waves from correlations of seismic noise and measuring their dispersion curves works rather well. However, further improving the accuracy of the noise based measurements needed to develop new high-resolution imaging and monitoring methods requires better understanding of the noise modal content and its evolution in space and time. Taking into account realistic distribution of noise sources is also necessary to be able to apply waveform inversion approaches to the noise correlations.

Seismic noise spectra contains two prominent peaks at 0.05 -0.1 and 0.1 -0.3 Hz called primary and secondary microseisms, respectively. The primary microseism originates from direct forcing of strong oceanic waves while the secondary microseism which is characterized by stronger amplitudes is produced at double frequency by a non-linear interaction of these waves as suggested by Longuet-Higgins [1950]. Both microseismic peaks are dominated by fundamental mode surface waves. It is currently debated whether the surface wave component of microseisms is generated primarily along coastlines (e.g., Friedrich et al. [1998]; [START_REF] Bromirski | The near-coastal microseism spectrum: Spatial and temporal wave climate relationships[END_REF]; Essen et al. [2003]; Schulte-Pelkum et al.

[2004]; [START_REF] Rhie | A study of the relation between ocean storms and the Earth's hum[END_REF]; [START_REF] Yang | Characteristics of ambient seismic noise as a source for surface wave tomography[END_REF]) or if it is also generated in deep-sea areas [START_REF] Cessaro | Sources of primary and secondary microseisms[END_REF]; [START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF]; Chevrot et al. [2007]; Kedar et al.

[2008]). At the same time, body waves were detected in the secondary microseismic band using dense seismic arrays (e.g., [START_REF] Backus | Wide-band extraction of mantle p waves from ambient noise[END_REF]; [START_REF] Toksöz | Microseisms: Mode Structure and Sources[END_REF]; Seriff et al.

[1965]; [START_REF] Iyer | vidence for the existence of locally-generated body waves in the short-period noise at the large aperture seismic array[END_REF]; [START_REF] Koper | Seasonal Anisotropy in Short-Period Seismic Noise Recorded in South Asia[END_REF]; Gerstoft et al. [2008]) and can

D R A F T December 1, 2009, 5:06pm D R A F T
be often associated with specific storms (e.g., Gerstoft et al. [2006]). Inhomogeneous distribution of noise sources is clearly revealed by the asymmetry of noise cross-correlations observed in both primary and secondary microseismic bands (e.g., [START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF]; [START_REF] Yang | Characteristics of ambient seismic noise as a source for surface wave tomography[END_REF]). According to Longuet-Higgins' theory, the generation of secondary microseisms is associated with the non-linear interaction of swells propagating in opposite directions. Such a configuration can be encountered in the coastal region where incident and reflected waves are likely present. This is the case that is considered as prominent by the seismologists based on the observations of the radiation by individual storms (e.g., [START_REF] Bromirski | The near-coastal microseism spectrum: Spatial and temporal wave climate relationships[END_REF]; Gerstoft and Tanimoto [2007]; [2009]). Although there is little doubt that a part of the ambient noise is related with the interaction of oceanic waves with the coast, it is not the only situation where waves propagating in opposite directions are encountered. [START_REF] Kedar | The origin of deep ocean microseisms in the North Atlantic Ocean[END_REF] used a wave action model to implement Longuet Higgins theory and found that particular regions in the deep oceans are potential sources of secondary microseism excitation. This is related to specific conditions of meteorological forcing associated with resonances of the water column.

Bromirski
Their results indicate that secondary microseisms can be generated in specific deep-water areas with one example in the Atlantic ocean south of Greenland.

To investigate the location of the sources of the background noise, we use seismological data averaged over long time series. The Rayleigh wave part of the noise in the secondary microseism period band rapidly attenuates with distance. It is therefore difficult to assess the locations of sources when the signal is dominated by the closest source, often the closest coast, that hides the remote sources. To overcome this difficulty, we use P wave at teleseismic distances recorded in continental environments.

D R A F T December 1, 2009, 5:06pm D R A F T
We compute noise cross-correlations for three seismic arrays located within continents in the Northern hemisphere. During summer months, when most of strong storms are located in the Southern hemisphere, the observed noise cross-correlations are dominated by arrivals at near-zero times. Polarization analysis clearly indicates that these arrivals are composed of teleseismic P-waves. We then use a beamforming analysis to determine precisely back-azimuths and slowness corresponding to these arrivals and to backproject them to the regions where the energy was generated based on ray-tracing in a global spherically symmetric Earth model.

Polarization analysis to detect teleseismic P-waves in noise crosscorrelations

In this section, we use the data of the Eastern Turkey Seismic Experiment (ETSE) that operated a temporary network of 20 broadband stations between October 1999 and August 2001 (Figure 1) [START_REF] Sandvol | Shear wave splitting in a young continent-continent collision: An example from eastern turkey[END_REF]).

Data pre-processing

Polarization analysis of noise cross-correlations requires preserving the amplitude ratio between components. Therefore, standard processing for the computation of noise crosscorrelations such as one-bit normalization and spectral whitening, can not be applied.

Instead, data are corrected for instrumental responses, resampled to 1 Hz and filtered between 0.01 and 0.3 Hz. A water level of 4 times the standard deviation of each record is used to decrease the amplitude effect of earthquakes on cross-correlations. Furthermore, for each station and each component, only daily records with mean energy smaller than the whole experiment mean energy are used for the polarization analysis. A10) for a plane P-wave propagating across a network of stations. Therefore, the eigenvectors are identical and the polarization analysis can be performed either on cross-correlations records at 2 stations or on 3-component record at a single station.

1.1.2. Noise cross-correlations D R A F T December 1, 2009, 5:06pm D R A F T
As suggested by [START_REF] Jurkevics | Polarization analysis of three-component array data[END_REF], the eigenvalues of the covariance matrix (λ 1 > λ 2 > λ 3 ) are used to compute the coefficient of rectilinearity R:

R = 1 - λ 2 + λ 3 2λ 1 . ( 1
)
which is equal to 1 for a rectilinear polarization. The eigenvector corresponding to the largest eigenvalue λ 1 gives the polarization angle and azimuth of the plane wave. The conversion from the polarization angle (ϕ) to the incidence angle (I) is obtained from the displacement equations for a reflected P wave at the free surface given by Aki and

Richards [1980].

We use a ratio between P and S waves velocities of : V P /V S = √ 3 for this conversion.

Polarization analysis of cross-correlations computed with ETSE data

We use only station pairs with distance larger than 50 km to compute the covariance matrix of the 3 cross-correlations ZE, ZN and ZZ. To prevent any influence of the Rayleigh wave (group velocity 3 km/s), we select time windows between -10s and 10s. Among the Using eigenvalues and eigenvectors of the covariance matrix, we compute the rectilinearity coefficient (equation 1), the azimuth and the incidence angle for every inter-station cross-correlation. The rectilinearity coefficient over the whole experiment is 0.84 ± 0.12, which shows that the polarization of the studied wave is almost linear.

Figure 4 shows particle motion for 2 daily cross-correlations and 2 station pairs. Particle motion is shown in the horizontal plane (ZN as a function of ZE) and in the vertical propagation plane defined by the measured azimuth angle (ZZ as a function of ZH) where ZH is obtained by the rotation of the ZE and ZN components of the correlations using the azimuth angle measured from the polarization analysis. The rotation can be computed after the correlation because no non-linear processing such as one-bit transform or spectral whitening has been applied to the data. The red dashed lines in Figure 4 display the azimuth and incidence angle obtained from the covariance method. We observe that the displacement is stronger on the vertical component than on the horizontal ones suggesting that the signal observed on the cross-correlations at near-zero times is composed of nearly vertically incident P waves.

We then estimate the incidence angle and the azimuth of the body wave detected from correlations of ambient noise records and investigate their possible seasonal variations.

Figure 5 shows the probability of occurrence of a given value of incidence angle (Figure 5a) and azimuth (Figure 5b) in a time period of 20 days evaluated from all daily records and station pairs. Figure 5a shows that we detect P waves with steep incidence angles during the whole experiment. It also documents a seasonal change of the incidence angle

D R A F T December 1, 2009, 5:06pm D R A F T X -11
from an average of 15 • in summer to 25 • in winter, with a more accurate measurement of the incidence angle in the summer than in the winter (larger probability of occurrence).

We observe the exact opposite in Figure 5b with better determined azimuths in winter than in summer, simply because the azimuth can not be evaluated for an almost vertically incident P wave.

The seasonal variation observed for the incidence angle is even clearer for the azimuth (Figure 5b). In summer, the average azimuth close to 0 [START_REF] Tanimoto | Seasonality in particle motion of microseisms[END_REF]). The precise location of the sources of the P wave component of the noise will be investigated in the following section.

Locating seismic noise sources with a beamforming analysis

To determine regions that generate these body waves, we perform a beamforming analysis of the noise cross-correlations using the whole network as an array. We use only vertical components where the body waves are mostly detected. When studying a single component, we do not need to preserve the amplitude and, for efficiency, pre-process the continuous data with spectral whitening and one-bit normalization to improve the signalto-noise ratio [START_REF] Larosse | Imaging from one-bit correlations of wideband diffuse fields[END_REF]). We analyze two frequency bands [0.05-0.1Hz] and [0.1-0.3Hz] corresponding to primary and secondary microseismic peaks, respectively.

D R A F T December 1, 2009, 5:06pm D R A F T

Beamforming analysis

Our time-shift beamforming analysis consists of decomposing the body-wave part of a wavefield recorded by a network into plane waves. If a plane wave defined by its slowness vector S reaches two stations A and B, the cross-correlation of signals recorded at these stations will be shifted by:

∆T AB ( S) = S • AB (2)
where AB is the vector connecting A and B. We approximate the network to be flat by

Locating P-wave noise sources on the Earth's surface.

In a next step, we project the results of the beamforming analysis on the Earth's surface.

Following the results of the polarization analysis of the near-zero-time arrivals at crosscorrelations and their fast apparent velocities, we assume that they are mostly composed of teleseismic P-waves. For a given slowness and back-azimuth we can back-project a D R A F T December 1, 2009, 5:06pm D R A F T seismic wave using a ray tracing in a spherically symmetric Earth's model. We suppose that diffracted and reflected phases are less energetic than direct and refracted waves and, therefore, the waves that we take into account are P, PP, PKP, PKiKP and PKIKP waves and we use the IASPEI91 tables [START_REF] Kennett | Traveltimes for Global Earthquake Location and Phase Identification[END_REF]) to relate the slowness with the source-receiver distance for all considered phases (Figure 7).

For every point on the Earth's surface, we identify all phases that may propagate from this point to the network location and determine their horizontal slowness. Therefore, for each position we determine the "energy" of all the waves that are considered from the horizontal slowness plane. The energy is evaluated from the function A determined by the beamforming analysis (equation 4). Finally, we select the phase corresponding to the maximum value of the beamforming map and attribute this amplitude to the projection at the considered geographical position. By repeating this procedure for all points on a 5 • longitude × 2.5 • latitude geographical grid we construct a map of what can be considered as probability density of noise sources during the considered period.

This process is illustrated for three geographical locations shown with stars in Figure 8b and ETSE Network (triangle in Figure 8b). White lines are the great circles and their projections in the horizontal slowness plane in Figure 8b Location of P-wave sources of the primary and secondary microseisms during different seasons are shown in Figures 10 and11. For every map we used correlations of one month of data. In the secondary microseismic band, regions that generate P-waves are well defined (Figure 10) and are mostly in deep oceans. Also, a clear seasonal migration can be observed with strongest sources located in the northern hemisphere during the northern hemisphere winter and in the southern hemisphere during the summer. In the primary microseismic band (Figure 11) we also observe a similar seasonal variation. However, uncertainties of the source location are much larger than in the secondary microseismic band. The reason for this is that the signal-noise ratio of the near-zero-time arrival in cross-correlations is much lower in the primary microseismic band than in the secondary D R A F T December 1, 2009, 5:06pm D R A F T microseismic band. Also, reliability of the array analysis degrades at longer wavelengths.

Despite these large uncertainties, it is clearly seen in Figures 10 and 11 that, in most cases, sources of primary microseisms do not coincide geographically with sources of secondary microseisms pointing to different physical mechanism of generation of these two microseismic peaks.

The maps of Figure 11 are similar with the results of [START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF] in terms of seasonnality. The location of noise sources based on the backpropagation of Rayleigh waves can not provide a resolution comparable to the one we achieved with body waves.

Nevertheless, even with P waves, the uncertainties of location of sources of the primary microseism do not allow to clearly conclude unambiguously that they are located in deep parts of the ocean, in near-coastal regions, or both. While strongest identified source areas tend to extend to deep oceanic parts, they also cover coastal areas.

Discussion and conclusion

Results of our polarization and beamforming analyses of the continuous noise records demonstrate that significant part of the microseismic noise is composed of P-waves generated by distant sources. These sources show a clear seasonality in correlation with the seasonal migration of the strong oceanic storms between the southern and the northern hemisphere suggesting that the observed teleseismic P-waves are generated by the interaction of waves produced by these storms with the seafloor. Location of P-wave sources in the primary and in the secondary microseismic bands do not coincide with each other indicating that these two peaks are generated in different regions and possibly by different physical process. P-waves are more easily identified in the secondary microseismic band than in the primary microseismic band. While we cannot exclude that this difference is

D R A F T December 1, 2009, 5:06pm D R A F T
related to the more efficient mechanism generating secondary microseismic P-waves than primary ones, a simple explanation of this observation can be related to difference in wave propagation. Strong noise sources generate both body wave and surface waves. For the latter, their attenuation is much stronger at higher frequencies. Therefore, surface waves in the secondary microseismic peak propagate much less efficiently over very long distances than in the low-frequency primary microseismic band. As a consequence, for distant noise sources, the relative part of the body waves in the recorded seismic noise is relatively high in the secondary microseismic band while the primary microseismic band remains largely dominated by surface waves, making observation of P-waves more difficult.

Using array-based processing of the teleseismic P-waves to locate regions generating strong microseisms has significant advantages relative to using surface waves (e.g., [START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF]). The latter yields only a determination of their back-azimuths at the array location while with body waves we can measure both backazimuths and slownesses that can be converted into distances. Therefore, we can locate the source regions more accurately with body waves than with surface waves.

We can compare our maps of the source density in the secondary microseismic band Overall our observations are consistent with the generation of microseisms by non-linear interaction of ocean waves propagating in opposite directions that create a pressure distribution on the seafloor at twice the frequency of the interfering waves [START_REF] D R A F T ; D R A F T Longuet-Higgins | A Theory of the Origin of Microseisms[END_REF]). Following [START_REF] Kedar | The origin of deep ocean microseisms in the North Atlantic Ocean[END_REF], this wave-wave interaction occurs in deep oceans and the geographical intensity of this interaction may be computed from oceanic wave action models. Moreover, the efficiency of the coupling between the interfering oceanic waves and the seafloor may depend on the depth of the water column, i.e., on the bathymetry.

As a result, an efficient transfer of energy from oceanic to seismic waves occurs over geographically very limited and specific areas. It is in particular interesting to note that the source area near Iceland seen in Figure 10 during October and January coincides with the strong source of Rayleigh wave microseisms computed by [START_REF] Kedar | The origin of deep ocean microseisms in the North Atlantic Ocean[END_REF] based on Longuet-Higgins's theory and oceanic wave action models.

These observations confirm that the source of secondary microseisms are not confined in the coastal areas as it is often accepted by seismologists. On average, the excitation of P waves by oceanic waves is stronger in the deep oceans. It does not mean, however, that there is no excitation along the coast, particularly when storms hit the shoreline.

In the future, locating sources of P-wave microseisms can be improved with using more networks better distributed over the globe. Observing and understanding these sources is important to validate predictive models of the seismic noise generation and distribution.

  D R A F T December 1, 2009, 5:06pm D R A F T vert travel times with different methods developed in context of earthquake or explosion based seismology.

Figures 2

 2 Figures 2 and 3 show cross-correlations between vertical noise records (ZZ) plotted with

  127,490 cross-correlation signals available for 671 days and 190 station pairs, 23,449 signals are selected based on signal-to-noise (see section 1.1.1) and minimum distance criterion. This rather small percentage (18%) is due to time-variable data availability and quality. D R A F T December 1, 2009, 5:06pm D R A F T Daily records including earthquakes or glitches are removed from the database due to our water-level amplitude filter.

  Figure6ashows results of the beamforming analysis applied to one-month cross-

  and 8a. For the location 1, the possible seismic phases are PP, PKP (branches ab and bc) and PKiKP ( red line in 8a) and the branch bc of the phase PKP corresponds to the strongest amplitude that is then selected for the projections. Similarly, for locations 2 and 3, possible phases are PKiKP and P or PP, respectively, and the latter correspond to stronger amplitudes and are used for the projection. The larger patch is projected as a P wave into the Indian ocean and D R A F T December 1, 2009, 5:06pm D R A F T as a PP wave into the southern Pacific. The smaller patch corresponds to a PKP wave originated in the vicinity of New Zealand. Maps of P-wave noise source densities corresponding to beamfoming results from Figure 6 are shown in Figure 9. Some source areas such as the region south of Africa during August 2000 and Northern Atlantic south of Iceland during February 2001 are well illuminated by both networks. This suggests that strongest sources of P-waves microseisms are seen by multiple networks distributed around the world. Therefore, we decided to combine observations from different networks to improve the accuracy of the location of main sources of P-wave microseisms. We used stations from three networks shown in Figure 1: 46 stations from Yellowstone park, 29 stations from ETSE, and 14 stations of the Kirgyz Seismic Network (KNET). We interpret the projection map obtained from every individual networks as a probability density and just multiply them to find the combined distribution.

(

  Figure 10) with results by Gerstoft et al. [2008] who applied beamforming to the noise records of the Southern California seismic network. We find similar source locations in southern Pacific and Indian oceans during summer months and in northern Pacific and Atlantic ocean during winter months. Using three networks simultaneously allows us to image source regions with higher reliability. One of most important consequences of this improved reliability is that we clearly see that strongest sources of P-wave microseisms are located in deep oceans far from coasts for the secondary microseism. Also, the observed D R A F T December 1, 2009, 5:06pm D R A F T source regions are significantly smaller than areas affected by significant wave heights.

Figure 1 .

 1 Figure 1. Maps of ETSE, Yellowstone and Kyrgyzstan networks

Figure 4 .SlownessFigure 6 .Figure 7 .Figure 8 .

 4678 Figure 4. Particle motion for 2 station pairs from noise correlations in the [-10s; +10s] time

  Figure 9. Projection results of the Yellowstone and ETSE networks for August 2000 and

Figure 10 .

 10 Figure 10. Seasonal variation of the location of P-wave seismic noise sources in the secondary

Figure 11 .

 11 Figure 11. Seasonal variation of the localization of seismic noise in the primary microseismic

  • shows that sources of the P waves are located south of the ETSE network. Winter noise sources are located north-west of the network as documented by azimuths close to 150 • . Those observations are consistent with seasonal changes in the behavior of seismic noise sources (e.g.[START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF];

Acknowledgments. The data of the "Eastern Turkey" and of the "Geodynamics of the Yellowstone Hotspot" seismic PASSCAL experiment as well as of the Kyrgyz Seismic Network used in this study were obtained through IRIS DMC. We thank Sharon Kedar and Frank Webb for helpful discussions and G. Moguilny for maintaining the Cohersis cluster. This work was supported by Agence Nationale de la Recherche (France) under contract ANR-06-CEXC-005 (COHERSIS) and by a EU project WHISPER.

Appendix A: Wave polarization analysis from cross-correlations

We consider a plane wave with amplitude A, wave vector k and pulse shape v(t). The azimuth (Az) is defined as the angle between the North and the projection of the wave vector in the horizontal plane. In the vertical plane, the incidence angle (I) is the angle between the vertical and the wave vector. Equation A1 gives the general form of a 3-component record of the signal that will be used for this demonstration :

where T is the arrival time of the signal at the station.

This plane wave is recorded at two different stations A and B at times T = t A and T = t B .

We compute three cross-correlations between the vertical component of station A and the 3 components of station B. They are :

where :

is the time correlation function and v(t) the complex conjugate of v(t).

The definition of the covariance between two signals (E and F ) in the time window [T 1 , T 2 ] is given by :

For North and East component records of station A, we select a time window including the

) and we compute the covariance defined by equation A4 to obtain :

where

We apply the same computation to all pairs of records at station A and we obtain the covariance matrix of the 3-component record :

where :

cos 2 (I) cos(I) sin(I) cos(Az) cos(I) sin(I) sin(Az) cos(I) sin(I) cos(Az) sin 2 (I) cos 2 (Az) sin 2 (I) cos(Az) sin(Az) cos(I) sin(I) sin(Az) sin 2 (I) cos(Az) sin(Az) sin 2 (I) sin 2 (Az)

We consider the definition of the covariance (Equation A4) and the definition of the crosscorrelations between 2 stations (Equation A2) to compute the covariance between ZE and ZN cross-correlations. We select a time window [T 1 ; T 2 ] of the correlation which includes the correlation of the initial signal :

where

The covariance matrix for cross-correlations is computed from Equations A2 and A4 :


Michael Eickenberg
email: michael.eickenber@nsup.org

Alexandre Gramfort

Bertrand Thirion

Scattering Transform Layer One Linearizes functional MRI activation in Visual Areas ICML Workshop on Statistics, Machine Learning and Neuroscience 2012

The scattering transform is a relatively new multi-layer and multi-scale signal transformation constructed to be invariant or tolerant to chosen transformations of the signal.

Here it is used in its simplest form to fit activation in visual areas related to the presentation of natural images. As it performs a spatial pooling operation over high-frequency phase invariant edge detectors, it mimics the spatio-temporal low-pass filter properties of the haemodynamic response oberved in functional Magnetic Resonance Imaging (fMRI).

Introduction

Only recently have biologically plausible models of primate vision become successful enough to be competitive with state of the art methods in the domain of object recognition. In order to address tractability problems, they are mostly inspired by Hubel and Wiesel's classic result [START_REF] Hubel | Receptive fields and functional architecture of monkey striate cortex[END_REF]) that the lowest level of processing is edge detection, the finding that visual cortex is hierarchically organized [START_REF] Felleman | Distributed hierarchical processing in the primate visual cortex[END_REF] and the fact that invariances and tolerances to certain transformations need to be built bit by bit along the hierarchy such that objects can be recognized despite differences in position, scale, orientation, illumination and pose. These ideas can be hardcoded into systems that implement the functionality of large groups of neurons and thus reduce the number of degrees of freedom left to be fitted to data if necessary. A common expression of this approach is the convolutional model, which alternatingly passes a linear filter (e.g. a linear edge detector, such as a Gabor filter) and a non-linearity over the input. In its 2007 version [START_REF] Serre | Robust object recognition with cortex-like mechanisms[END_REF], the HMAX model achieves state of the art object recognition performance and humanlike performance and errors in rapid object recognition tasks. It is implemented as an alternation of linear operations and pooling over space and scale using a maximum operation. The convolutional neural networks in [START_REF] Lecun | Convolutional networks for images, speech, and time-series[END_REF] achieved state of the art performance on the recognition of hand-written digits.

Compared to electrophysiology, fMRI is a relatively new method to localize brain function. Albeit orders of magnitude coarser in spatio-temporal resolution and indirect in the nature of the signal, it permits the acquisition of the activity pattern of all the visual (and other) regions of a brain while a subject is viewing an image. Recently, fMRI brain activity of subjects viewing natural images has been shown to provide enough information to be able to identify the stimulus among a very large set of images [START_REF] Kay | Identifying natural images from human brain activity[END_REF]) using an encoding model, and even perform a reconstruction of the stimulus using a "bag-of-images" prior and Bayes' theorem [START_REF] Naselaris | Bayesian reconstruction of natural images from human brain activity[END_REF].

The advances in convolutional models for object recognition have led to an increased interest in their mathematical properties. Some degree of invariance is one of them and constitutes a key property of the signal processing in the brain. The scattering transform [START_REF] Bruna | Classification with scattering operators[END_REF][START_REF] Mallat | Group invariant scattering. to appear[END_REF] is a recent development and is built to create invariance to any chosen transformation group in a mathematically rigorous way. It consists of an alternation of an analytic directional wavelet transform and complex modulus followed by smoothing. The complex modulus is a non-linear component which adds stability by losing the phase of the transformation and prevents the smoothing from losing information (by averaging phase dependent edge responses). It turns out that the local averaging properties of the scattering transform may be well reflected by fMRI voxels.

In this paper, the scattering transform will be briefly introduced and the analogy between its first layer and visual fMRI voxel properties will be established. The first scattering layer will be compared with the encoding model of [START_REF] Kay | Identifying natural images from human brain activity[END_REF], which consists of Gabor filter energies, and contrasted in predictive power by cross-validation.

Notation

For u ∈ C the complex modulus is noted |u|. Real signals are noted f and wavelets are ψ with scaling function φ.

The scattering transform

The scattering transform, introduced in [START_REF] Mallat | Group invariant scattering. to appear[END_REF], is a signal transformation based on analytic wavelets which is constructed with the goal to establish a representation invariant or tolerant to any chosen transformation group on the signal while keeping its essential characteristics. By default, in images, the transformation group is translations and the essential characteristics are texture or object features. This property is enforced using a cascade of directional complex wavelet transforms followed by a complex modulus. At each stage, a smoothed (low-pass filtered) version of the intermediary result is sent to output, while the high frequencies are kept by the subsequent wavelet transform. The complex modulus acts as a convolution in Fourier space, which picks up correlations at a given frequency difference, thus bringing high frequency interactions within the Fourier support of the wavelet down to their frequency difference on each layer.

Here smoothing will be implemented as a local spatial average, building local translation tolerance.

Let Γ be a set of angles and ψ γ , γ ∈ Γ a directional analytic wavelet, e.g. a Gabor filter. With φ J as the scaling function, we scale for j = 0 to j = J as follows:

ψ γ,j (x) = 1 2 j ψ γ (2 -j x).
The first layer of the scattering transform of signal f is then

W γ,j f = |ψ γ,j * f |
The values of the first layer are smoothed

S γ,j f = |ψ γ,j * f | * φ J
and are output as such or globally averaged. The sec-ond layer and its output are as follows

W γ1,j1,γ2,j2 f = ||ψ γ1,j1 * f | * ψ γ2,j2 | S γ1,j1,γ2,j2 f = ||ψ γ1,j1 * f | * ψ γ2,j2 | * φ J
This can be continued and is further explained in [START_REF] Mallat | Group invariant scattering. to appear[END_REF]. The signal energy does not diverge and for appropriate wavelets the transformation is unitary. In analogy to the outputs of layers 1 and 2, the 0th layer output reads as follows:

S 0 = f * φ J

Scattering layer 1 and fMRI voxels

The early visual area V1 contains simple and complex cells that perform phase dependent and phase invariant edge detection. A contour can be defined by its location in the visual field, its orientation and its spatial frequency. A large number of these detectors is arranged in a retinotopic manner across V1. Retinotopy preserves local neighbourhoods, thus neurons that are close to each other will receive visual information from nearby positions in the visual field. This means that one fMRI voxel in V1 contains many neurons tuned to the same region in the visual field, while orientations and scales can be different. The hemodynamic response measured with fMRI can be seen as a spatiotemporal low pass filter. Since all the neurons inside the voxel contribute to its response, we can attempt to model it by locally averaging complex wavelet moduli of all scales. This method keeps high frequency energy, but blurs it out to the scale of the voxel. In fact, it can be viewed the first layer of the scattering transform. A similar principle is described in the DAISY transform (Tola et al., 2007) -a method for extracting scale invariant features that differs from the classic SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] in that it takes local directional derivative averages instead of histograms of derivative orientation.

Methods

We use Morlet filters, i.e. Gabor filters with 0 DC offset in the cos part. This is achieved by choosing C > 0 such that the integral over exp -

1 2 x T Σ -1 x (cos(k T x) -C)
is zero. We use these filters throughout the experiment in order to permit a comparison between the two transforms.

We compared the performance of the first layer of the scattering transform on the [START_REF] Kay | Identifying natural images from human brain activity[END_REF] dataset with the original filter pyramid that was employed. To this end, we implemented the pyramid according to the information provided in the supplementary material to the article: At 8 orientations and 5 scales, cos and sin type real, isotropic Gabor-like Morlet filters are placed at discrete locations on the stimulus image as follows:

The middle of the image for the largest scale, at a spatial frequency of one cycle per image size. For the next scale the middles of the four sub-images of half the side length obtained by cutting through the middle of the image are occupied by filters one octave smaller, at two cycles per image size. The next scale is obtained by cutting each of the previous sub-images into 4 subsub-images and placing a filter in their center points. This process is continued twice more in order to end up with 256 filter locations on the smallest scale. The cos and sin parts of the filter response are combined by using the square root of the sum of squares, leading to a complex wavelet modulus as the transformation output. Adjustment parameters were the number of scales, J, chosen from {3, 4, 5, 6}, and the size of the first wavelet, determining the size of all the wavelets in the pyramid, chosen between a scale factor λ s ∈ { 3 4 , 1, 5 4 } times the original size. Using the same type of wavelet, the scattering transform was performed using J ∈ {3, 4, 5, 6} scales and 8 orientations. The output smoothing is defined by the largest scale through the convolution with φ J . Hence, the size of this low-pass filter must be in accordance with the size of the receptive fields of the voxels, and especially not too large.

The dataset of [START_REF] Kay | Identifying natural images from human brain activity[END_REF] consists of 1750 training stimulus images and 120 test stimulus images of size 128x128 pixels along with preprocessed fMRI voxel responses for two subjects of around 25000 voxels for each image. Minimally preprocessed fMRI data is also provided, but here we use only the preprocessed responses from both subjects.

In order to evaluate the linear predictive performance of the transformations, we perform a ridge regression on the transformation output as the data and the voxel response as the target value. We use 5-fold nested 5fold cross-validation on the training set and predictive r 2 as the scoring function. The test set is not used for evaluation, as it is too small for our purposes and bears the risk of overfitting when used for parameter setting.

Results

For the Gabor pyramid the optimal settings measured by predictive r 2 are J = 5 scales and scale factor λ s = Figure 1. Scatterplots comparing predictive r 2 scores on subject 1. Each cross represents a voxel. The score using the original filter pyramid is on the x-axis. The y-axes are the scores using the scattering transform with different maximal scales J. Top left: J = 3, top right: J = 4, bottom left: J = 5, bottom right: J = 6. At higher scales the output is too coarse to fit voxel activity well. Scale J = 4 has the best results. Scale J = 3 has slightly worse results, possibly due to a higher number of coefficients in the regression.

1, as used in the original paper. A 6th scale adds 8192 coefficients to the 2728 of the setting J = 5, thus adding a lot of data which does not seem to be relevant to voxel prediction. On the other hand, J = 4 and less provides too few and too coarsely resolved coefficients. Hence we restrict ourselves to comparing the results obtained with these optimal settings to those of the scattering transform.

All figures show data from subject 1, but the presented results equally hold true for subject 2. In Figure 1 we provide scatterplots comparing the linear predictive performance of the original model (x-axis) to the scattering transform first layer (y-axis) at different numbers of scales. For 3 and 4 scales the point cloud clearly lies above the diagonal. A Wilcoxon signed rank test confirms this with p ≈ 10 -120 and p ≈ 10 -200 respectively. The scattering transforms using 5 and 6 scales are nevertheless comparable to the original filter pyramid. In these cases, the output smoothing mechanism yields a response image size of 8x8, respectively 4x4 pixels per scale and orientation, which is too coarse to be able to resolve the receptive fields of typical voxels. As few as 3 scales in Fig. 1 top left are sufficient for the scattering transform to fit the voxel responses very well. In this case the output images are smoothed and downsampled to 16x16 pixels and are in accordance with the size of voxel receptive fields. This shows that the voxel activity can be fitted to edge filter responses of very fine scale, pooled over the region of its receptive field.

Figure 2 shows the receptive fields of two selected voxels for subject 1. They are visualized using different techniques: On the left there are the receptive fields due to the scattering transform, generated by taking the mean of the response over all scales and orientations. On the right, in order to visualize a receptive field using the original filter pyramid, we multiplied each resulting coefficient by the envelope of the corresponding filter. The discrete spatial nature of the pyramid becomes clearly visible. Both methods localize the receptive field in the same area and the signal energy is equally distributed over all orientations (not shown).

Discussion

In this paper we introduced the scattering transform as a candidate to model visual voxel activity. We compared it to the original Gabor pyramid and showed that fitting the size of a pooling region to the size of a voxel's receptive field yields higher predictive r 2 scores than fitting the size of the filter to the size of the voxel's receptive field. This gives rise to the hope that the subvoxel activity, which is composed of many visual neurons with receptive fields smaller than the one of the voxel, can be captured by attempting to model them and then pooling the result onto voxel size. Evidently, a rise in predictive score cannot be taken as evidence of having found a better model, but nevertheless, further inquiry into modeling of subvoxel activity followed by pooling as a means to model fMRI activity patterns seems called for.

Figure 2 .

 2 Figure 2. Receptive fields of two selected voxels using J = 4 scattering on the left and the original filter pyramid on the right. Receptive field size is on the order of 1 8 of the side length. At least 1 16 of the side length is needed to resolve it well.

  Tola, Engin, Lepetit, Vincent, and Fua, Pascal. A fast local descriptor for dense matching. Technical Report École Polytechnique Fédérale de Lausanne, CVLAB, 2007.

Acknowledgments