Arnadou Ndtaye

Francis Fleurat-Lessard

Diagnosis and grading of wheat grain initial quality by a computerised decision support system

This paper descnbes the qualitative reasomng used to assess the imtial quahty of gram The mam problem was the representation of a heterogeneous knowledge m the same homogeneous calculus space, and the development of reasonmg methods to handle a formal knowledge base on wheat, maltmg barley, maize or any other cereal gram. This modifiabihty implies a user mterface with a separation of the static descnption of the mterface (corresponding to the knowledge base) from the mterface engme which mterprets the static descnption A prototype of knowledge based system, QuahS, was built up on Wmdows platform.

Introduction

In these last ten years, there have been a great mterest on usmg the knowledge based systems to assist the management of stored gram. Several pest control systems have been developed to support the control of stored gram pests (msects and mites) (Flmn et al., 1990), (Wl1km et al , 1991), [START_REF] Longstaff | A deCision support system for pest management m the Australian central gram-handling system[END_REF], [START_REF] Pasqual | Development of a Prototype System for Identification and Control of Insect Pests[END_REF], [START_REF] Jones | Development of an expert system for pest control in tropIcal grain stores[END_REF]. Today, current research focuses on preservation and gradmg of gram mitral quahty [START_REF] Ndiaye | Research on an expert system for appropriate management of the quahty of stored gram for food and feed processing[END_REF].

Gram quahty preservation has traditionally been performed by store-keepers who rely on measurements and observations on gram and its millmg products, and on imphcrt knowledge gamed through scientifrc results, common sense and job expenence. Each store-keeper uses his own method to evaluate and grade the initial quahty. Gram quality preservation implies both a correct assessment of gram uutial quahty and of its fmal use The mam difficulty encountered m rapid gradmg hes m the accurate evaluation of the gram mrtial quahty (Wngley et al., 1994), [START_REF] Maier | Gram Quality Task Force. Quality Grain Needs TLC[END_REF] One of the most commonly used method is gradmg accordmg to gram vanety (Morns & [START_REF] Morris | A computer-aided approach to the evaluation of wheat grain and flour quahty[END_REF]. But it can happen that the properties of 1 INRA, Laboratoire des Insectes des Denr's Stock's, BP 81, F 33883 Villenave d'Ornon, France EmaIl: amadou ndiaye @ bordeaux mra fr a vanety of wheat may change from one farming area to another and from one year of harvest to another We have defmed a qualitative model to assess grain wheat nutial quahty based on the measurements and observations of the gram standard charactenst~cs In this paper, we present the quahtatrve representation and reasonmg on the heterogeneous knowledge on gram preservation dunng storage. The knowledge base system QualtS is also descnbed.

The initial quality of grain

The mitral quahty of gram is a new an complex concept, there is no model for its quantification In the literature, models of quahty changes dunng storage relatmg the maximum safe storage penod or loss of dry matter to the temperature and moisture content or carbon dioxide production are proposed The mam charactenstics of the gram quality studied were the loss of germmative Vigour [START_REF] Armitage | Pest control by coohng and ambIent air drying. m Flannigan B. SpoIlage and mycotoxius of cereals and other stored products[END_REF]Favier and Woods, 1993), moulds growth [START_REF] Fraser | Airflow requirements for drying gram with ambient and solar-heated aIr in Canada[END_REF][START_REF] Armitage | Pest control by coohng and ambIent air drying. m Flannigan B. SpoIlage and mycotoxius of cereals and other stored products[END_REF][START_REF] Latif | Respiration model for heat and moisture release during grain storage Amencan Society of Agncultural Engmeers[END_REF], and msect pests dynamics which are specific qualities charactenstics. The germinative Vigour is a technological characteristic, and the mould growth and the msect pests dynanucs are samtary & safety charactenstics. Neither all the characteristics of the quahty of gram nor their interactions have been modelled: for example the mfluence of the protem content on the gram quality is not modelled

The quality of gram is generally defmed as the suitability for the end-use of the gram: "the word "quahty ' means surtabihty for the specific process or utrhsation for which it'S destined" (Wngley et al , 1994). To estabhsh a sound qualitative or quantitative model of the quahty of gram, it is necessary to distmguish the quahty to the end-use Each utilisation implies a given quahty of gram. This required quahty must be available at the time of utilisation to make it possible. Therefore, to preserve the stored gram quahty compared to its possible end-uses, it is imperative to determme These vanables are heterogeneous considering the umts of measurements (seconds, kh -Ilhl, %, etc.) and the types of values (numerical or lmgurstic). They were interrelated by cause-effects relations.

The causal relations between the variables

The qualitative reasoning to assess the grain initial quality

The quahtatrve reasonmg is based on the theory of Quahtative PhYSICS whose about 20 years old field of Artificial Intelhgence (Hayes, 1979) The Quahtative

PhYSICSenable us to cope WIth mcompletely known dynamic system usmg heterogeneous data and the reasomng on these data does not necessary require the use of quantitative 2. to integrate in the same computer system gram characteristics (measurements and observations on gram and ItS nulhng products), and knowledge from scientific results, common sense and Job experience. 3. to build up a modifiable knowledge base system to be used in different environments: a silo where the nutritional charactensucs are not taken mto account, and a silo where they are; or a silo where they have an mstrument to measure the wheat bakmg quality and another where they have not etc ...

Qualitative modelling

A static causal graph has been developed as the structural model of the causal relationships between the vanables It was an oriented graph where the nodes were the vanables-: the measurements, the observations and the calculations -, and the onented lmks (arcs) were their relationships (Figure 1).

The qualitative representatives of measurements and observations, and the calculations were defmed in the same discreet quahties space Work done on qualitative reasomng, and particularly the Dual formahsm developed by [START_REF] Guerrm | Duahstic algebra for qualitative analYSIS[END_REF], seemed particularly efficient to represent the qualities space Applying this formahsm, it became possible:

• to define the qualitative representatives of measurements and observations, and the calculations mto the same quahties space . to determme the computmg operators on tlus quahties space . to defme a vocabulary to descnbe the gram conditions. The qual itiee space Each vanable was defmed mto a specified validity domam Here, for every vanable of any type, ItS vahdity domain could be divided mto a maximum of seven mtervals of values, representmg the relevant qualitative differences For example the moisture content (me)

The compuiuut operators

The operators which allow combmation of human-experts knowledge m the quahties spaces were defmed. These operators were built up as ad hoc decision tables. They enabled us to put mto words the projection laws for measurements and observations into the quality spaces, and the combmation laws withm these spaces. These laws were aSSOCla ted with weighmg factors ranging from zero to one (Tables 3 and4).

Grain too humid

Grain heated

Grain moderately heated

In spite of its simple graphical representation, the causal graph is a powerful tool for modelling a complex and modifiable system. The stored grain ecosystem is a complex system where several variables interact and influence the grain quality.

In this domain, it happens that new instruments are developed in order to measure new characteristics improving decision making for grain quality diagnosis. In this case it would be very easy to add a new node in the causal graph: and if one test is not longer used, it will be easy to retract a node from the graph.

Qualitative reasoning

We established the abstract operations to describe functionally the causal relations between the variables, a method to propagate the weighing factors from the human experts in the causal graph determined, and a nonmonotonic reasoning to handle the modification of the value of measurements and observations.

The abstract operatiems

The determination of abstract operations was based on the distinctive features of meaningful groups and the nodes and links of the causal graph.

MC SW vi I m h vh wh I 1 1 0.8 0.9

1 1 m 1 1 1 0.8 1 0.9 h 1 1 1 1 1 1
Meaningful groups: we distinguish the meaningful groups which only give information on other meaningful groups Nodes and links: we distinguish the nodes without predecessor (leave nodes) to the nodes with predecessor.

(
The types of nodes are: measurement (mes-X}, observation (obs-X), and calculation X. There are two types of links:

influences and gives-information-on.

The criteria of our reasoning on a node were the number and the types of its predecessors and the types of links.

From these criteria and the distinctive features we defined six abstract operations to calculate the qualitative value of (the inputs), a node to be calculated (the output). and an operator to be executed in order to get the value of the node to be calculated.

• a Projection was a one argument operation which declares a projection of a measurement or observation in the qualities space; its mput was a leave node, and its output was a calculation. • an Induction was a one argument operation which declared an mduction of a calculation node from another calculation node. The relation between the prohferation of microorganisms (PmO) and the ergosterol content (EC) in the Samtary & safety characteristics meaningful group (figure 1 (5» illustrate an Induction: EC which is mfluenced by PmO is measured to have information on PmO • a Reduction was a two argument operation which declared a reduction of two calculation nodes mto another calculation node. For example, the baking quality (BQ) is a reduction the baking index (W) and the extensibility mdex of bread dough (P IL) . • a Prediction was a two argument operation which declared a prediction of a calculation node from two other nodes.

The risk of proliferation of nucro-orgamsms (PmO) is predicted from the temperature of gram (q) and the relative humidity of air (rh). • a Deduction was a one or two argument operation which declared a deduction of a calculation node from one or two other calculation nodes. The hidden form of msects (HF) could be deduced from the visible living form of insects (LF) and the visible dead form of insects (DF) • a Critic was a one or two argument operation which declared a improbability of the value of one calculation node at a given moment or an inconsistency between two DF VISible dead forms ofmsects EC Ergosterol content ,....--------------"""1-HF Hidden forms of insects The operations were used to represent m the knowledge base the causal relations between the vanables (Figure 2) and to indicate how to use the computmg operators. The operations were meta-rules which indicated how to generate new knowledge from raw data and results: they were used to represent m the knowledge base the causal relations between the vanables and to indicate when to use the computing operators (Figures 2 and3) .

Propagatwn of uncertauitu

In tlus situation of an expert domain, events were not necessarily repeatable, and the exact mathematical form for the uncertamty and the relative frequency of errors were unknown. The human expert had given weight factors to his decision rules in associatmg them With a subjective probability which came from his own expertise. The weighing factors associated within the computmg operators were propagated through the calculus space to determme the certamty factor of each calculation node [START_REF] Ndiaye | Diagnosis and Grading of Grain Initial Quality[END_REF] using the followmg formula: CF=CxWF with CF the certamty factor to compute, WF the weighing factor of the operator to be executed, C the confidence level lmked to the likelihood of an event, symbolised by the presence of one or several objects m the system status space: A nonmonotonic reosonmq A nonmonotonic reasonmg was necessary to allow the user to SImulate the decrease or mcrease of the measurement values (temperature, moisture content, ergosterol content, etc. ). Every modification of the value of a measurement modifies the value of every node directly or indirectly mfluenced; and consequently modtfies the assessment of the quality of gram [START_REF] Ndiaye | Preservation of stored grain quality, diagnosis and grading of wheat gram initial quality[END_REF].

QualiS: a prototype of a modifiable knowledge base system

QuahS IS an on-going expert system project for the preservation of gram initial quality durmg the storage penod QualiS was structured m three mam parts: a modifiable user mterface, a knowledge base and an mference engine. ThIS structure allows us to build up a modifiable knowledge base system to be used m different silo environments.

The graphical user mterface allowed the user to Ignore totally the reasonmg engme language, and to exchange mtuitrvely WIth It [START_REF] Pon | Relisation de l'mterface graphique modifiable d'un system expert de diagnostic et de prervation de la qualite du grain[END_REF]. The knowledge base contained the facts and rules representmg the human experts knowledge: the variables of the gram qualIty domam (nodes of the causal graph), qualitIes spaces, operators, weIghmg factors and defined vocabulary are facts; the operations, the models to propagate the weighing factors are rules. The mference engme integrated the formal procedures and functIOns to compute the knowledge.

The modifiable graphical user interface

The graphical user mterface has been Implemented in 1930 Allegro Common LISP [START_REF] Steele | Jr 1990 Common LISP the Language, 2d edition[END_REF]) For easy learning and usmg by a non computer specialist, this interface used Wmdows standard controls, as buttons, lists, etc. (FIgure 4) . These controls were organised followmg the classical style of Wmdows applications. This mterface could be completely controlled WIth the mouse, and many ergonomic details made It easier to use: as spin-boxes to help entenng numbers, help tool tips indicatmg the meaning of every graphical object, etc.

As the system may be changed, to follow the progress m measunng tools and techniques, or to work WIth another type of cereal gram, the graphical user mterface had to be built up in order to facilitate the software settmg up. The descnption for the wmdows and data was loaded at run time.

ThIS descnption was divided in two parts: the static part (the descnption in a specific language of the windows and their contained objects) and the dynamic part (the dynamic relations between the objects).

The static part: descnbed all the objects mvolved: the displayable graphical objects (the windows, buttons as OK button for example, the spm-boxes to assist the measurements entry, spm-check-boxes'' for the data entry of non necessanly available measurements as the ergosterol content, etc) and the static relations between graphical This part was wntten in a language to be used by non computer specialists, It was very simple, yet complete. The following example presents a language close to French (figure 5). In fact, this language was translated at run-time into English, and internally treated in English. This on-thefly translation was achieved by the use of a dictionary file which provides, for each keyword, a conversion between English and the user language.

Each action onto the graphical interface was associated to a fact, specified in the description. At run time, these facts were transmitted to the inference engine, which became informed of what the user did.

The dynamic part: described the dynamic relations between the graphical objects. This description allowed actions to start when some conditions became true. The knowledge base contained a declaration of all the existing dynamic relations between the user interface objects. These relations were represented as couples (conditions, action); so when the conditions of a dynamic relation became true, an order was given to the user interface to execute the corresponding action. For example when the number of visible living or dead form insects was more than zero, the push-button (RartitIon des insectes ) associated to the Identification of insects window was made enable; or when the cancel button of a unndou: was pushed, all the facts coming from the objects contained in this window were retracted, the facts derived from them to and the window hidden.

The knowledge base

The knowledge base has been implemented in CLIPS (Giarrantano and Riley, 1994). It was a set of declarative knowledge on preservation of stored grain quality, actually the baking wheat, and the knowledge on the dynamic part of the user interface [START_REF] Ndiaye | Preservation of stored grain quality, diagnosis and grading of wheat gram initial quality[END_REF].

The inference engine

The inference engine has been Implemented in CLIPS. It computed every knowledge base which was built up using the specific representation described in [START_REF] Ndiaye | Diagnosis and Grading of Grain Initial Quality[END_REF], and monitored the display of the main windows of the user interface [START_REF] Ndiaye | Preservation of stored grain quality, diagnosis and grading of wheat gram initial quality[END_REF].

Discussion

Evaluating the system QuahS qualitative reasoning approach represents a significant advance in existing computer knowledge based systems for preservation of stored grain. The departure from pest control to quality preservation is a big change for the stored gram conservation This is clearly an advance for the grain storage industry, which previously could deal with the total quahty of grain only m very specific cases. QualIS allows store-keepers to assess the nutial quality and grade of gram at delivery. It further provides the grain conditions, the complementary tests to be done to refme the grain quality, the incompatible uses of the grain, the grain quality altera tion risks and the storage operations to be done to preserve the gram quality (Figure 6).

1) «Identificateur

On going work

In this version of QualiS we are developing a results display module of the user interface allowing the user to easily ask questions on the results. When the results are displayed in natural language as 'The grain quality is considered as very good with a suggestive evidence of 0.95' , where the keywords are underhned, the user would be able to ask questions on these keywords. For example What is? for the keyword grain quality, or Why? for very good. This will be done with a simple mouse chck on the keyword; and then a floating menu will be displayed, allowing the user to choose, m the list of available questions, the question he asks.

In the future version of QualiS, we will develop a planning and a monitoring modules as parts of the reasoning engine; also we will develop an interface for the maintenance of the knowledge base by the human experts.

The goal of the planning module is to generate an optimal technical route for the storage of grain; and the monitoring module will predict and control the grain quality evolution in time. If it is a significant gap between the prediction and the grain condition, the storage technical route will be re-planned after a diagnosis of the causes of the gap. -II Y a une forte presence d'insectes sous forme libre morts (0.9) II est consellle de faire Ie test complementalre de recherche des formes cachees d'insectes (0.9) II est recommande de lui Appliquer de I'lnsecticide remanent (0.9) -la temperature du grain est trop chaude (0,9) II est recommande de Ie Refroidir (0.9) -Ie grain est echaude mediocre (0.85) II est recommande de ne pas I'utlllser en melange avec un autre gram (0,85) II est recommande de Ie Nettoyer (0.85) -Ie grain est contamine par les mycotoxines (0.9) II est recommande de ne pas I'utlllser i I'exportatlon (09) II est recommande de ne pas I'utiliser en consommation humalne (0.9) II est recommande de ne pas I'utiliser en amidonnerie (0.9) II est recommande de Ie declasser en alimentation ammale s'll n'est pas toxique (0.9) II est consellle de faire le test complementaire de recherche de genres et d'especes de micro-oroanismes ro.91

  before the storage what is the gram irutial quahty and grade m function With the gram charactenstics The baking wheat characteristics Store-keepers are used to evaluatmg the gram quahty through measurements and observations on gram and Its milling products. Measurements are numenc values from standard tests (temperature, moisture content, protem content, ) and observations are Imguistic values from sensonal tests (colour and smell) or are given by the gram producer (variety, area of production and year of the harvest). We have distinguished seven meanmgful vanables groups and determined five relevant calculated variables used to represent the four main specific qualities and the global quality of grain: • Grain identity variables group: the variety of the gram, its area of production and the year of harvest are given as linguistic values by the producer to the store-keeper. • Gram aspect vanables group: the colour and smell of grams are sensorial tests done by the store-keeper, the results of which are lmguistic values. These tests could be done automatically by Image analysis. arr and temperature of gram, and the carbon dioxide release rate in the gram environment are more or less easily measured by the store-keeper, the results of which are numenc values. • Intrmsic & physico-chemical characteristics vanables group: the moisture content, specific weight and impurity rate of grains, and the weight of 1000 grams are measured by the store-keeper, the results of which are numeric values • Samtary & safety characteristics vanables group: the presence of insects (VIsible livmg and dead forms, and hidden forms) or mites and the proliferation of microorganisms in gram, the mycotoxins and the pesticide residues on the gram, and the heavy metals and radioactivity m the grain can be assessed by the storekeeper or a speciahsed laboratory, the results of wluch are numeric values. • Technological charactenstics vanables group: the protein content, baking quality of wheat and rye flours (baking index Wand extensibility index of bread dough P IL), enzymatic actrvity, quality of the starch, germmation vigour and gram hardness can be measured by the storekeeper or a nuller, the results of which are numeric values. • Nutritional charactenstics vanables group: the quality of lipids, digestibihty of protems, and nutritional value can be measured using speciahsed mstruments by the storekeeper or a feed plant laboratory, the results of WhICh are numeric values. • Global quality charactenstics vanables group: the mtnnsic & physico-chemical quality, sanitary & safety quality, technological quality and nutritional quality are the specific qualities, and the gram quality can be associated by the store-keeper as qualitative values.

  problems we mtended to solve were: 1 to mtegrate in the same calculus space all the characteristics of gram, the pertment vanablcs. ThIS problem mcludes the choice of an appropnate 1925 vocabulary to describe the gram conditions.

Fig. 1 .

 1 Fig. 1. Causal graphs of variables: (0) represents the causal relations of the meaningful groups, and the (1) to (8) represent the causal relations of the vanables of each mearungful group

  Grain aspect and Grain identity), the meaningful groups which influence the future of the grain quality (Environment factors), and the meaningful groups concerning the grain conditions or characteristics (Intrinsic & physico-chemical conditions, Sanitary & safety characteristics, Technological characteristics, Nutritional characteristics and Global quality characteristics) : • The variables included in the groups Grain aspect and Grain identity were only used to supply the lack example, if we could not measure the ergosterol content to induce the qualitative proliferation of micro-organisms, the colour and smell of grain should be used to deduce less precisely the qualitative proliferation of micro-organisms. , The variables included in the group Environment factors were used to predict the risk of grain damaged by the micro

  (Deduction. Prediction and Critic). An abstract operation was a relation between one or two nodes

Fig. 2 .

 2 Fig. 2. A srmplified mmal quality assessment process for bakmg wheat

Fig. 4 .

 4 Fig. 4. The mam wmdow of the graphical user mterface

  Fig.5. Static part: 1) descnption of spin-box temperature; 2) description of the mam wmdow which contains group-boxes (etat-phYSIC0-chlIIl1que qualrte-technologique etat-sarutaire ) and a group of push buttons (groupe-des-boutons-fenetrepnnopale), Every group-box contams the descnption of graphical objects, and the group of push buttons contains the descnpnon of push buttons.

  et securlte Mauvals(e) (0 7of) _____ the mcompatible uses of the gram and/or the L'etat du grain est Ie suivant .

11 1+Fig. 6 .

 116 Fig.6. The results are presented through four parts: the gram quahty. the grain conditiOns and its lmked adVice. the grade of the gram and the remarks (degradation nsks and mconslstencles)

Table 1 .

 1 List of the relevant vanables m the space of the quality of bakmg wheat.

	relations.	The	qualitative	approach	completes	the	qualitative knowledge, but does not mtend to replace it in a
	quantitative one by allowmg the modelhng of human experts	well-known numerical problem space
	~eaningful group		Vanable				Abbrevia tion	Value from
	Gram Identity		variety				Va	Observa tion
				area-of-production			AP	Observa tion
				year-of -the-harvest			YH	Observa tion
	Grain aspect		colour				C	Observation
				smell				S	Observa tion
	Environment factors	relative-humidrty			RH	Measurement
				tempera ture			El	Measurement
				carbon-dioxide-ra te			C~Measurement
	Intnnsic &			moisture-content			MC	Measurement
	physico-chemical		specific-weight			SW	Measurement
	characteristics		impurity-rate			ImpR	Measurement
				weight-of -1000-grams			WlOOO	~easurement
	Samtary & safety		visible-hvmg-form-of-msects		LF	~easurement
	charactenstics		visible-dead-form-of-msects		DF	Measurement
				hidden-form-of-insects			HF	Measurement
				presence-of-insects			PIns	Calcula bon
				ergosterol-content			EC	Measurement
				prohfera non-of -micro-orgamsms	PmO	Calcula tion
				mycotoxms-			MT	Measurement
				pes ticide-residues			PR	Measurement
				presence-of -mi tes			P~It	Measurement
				heavy-metals

  ISdefmed between 8 and 22 %: If me ~12% the gram is very very dry, If 12 < me ~13% the gram is very dry, if 13 < me ~14% the gram ISdry, If 14 <mc ~15% the gram ISnormal, If 15 < me ~16% the gram IS humid and If 16% < me the gram IS very humid and could not be preserved m this condition.The qualities space QS composed by seven symbohc elements 1 vvl, ol , 1, m , h , vh, vvh I were used to represent the significant mtervals of values of each vanable (Table2). The elements were orderedvvi < vI < I < m < h < vh < v'vh -and elementary operators were defmed(Guernn, 1995) The qualities space allows us to mampulate m the same space vanables With very different interpretations and semantics; this property IS very Important for our purpose

							Global quality
							charactensttcs
					Nutntional
					characteristics
			6		7		8
	B	Temperature	ImpR	Impunty rate	Qn	Nutntional quality
	AP Area of production	LF	Visible Irving form of msects	Qss	Samtary & safety quality
	BQ baking quality	MC	MOIsture Content	Qt	Technological quality
	C	Colour of grams	MT	Mycotoxms	QL	QualIty of lipids
	C~Carbon dIOXIde rate	NV	Nutntional value	QS	QualIty of the starch
	DF Visible dead form of insects	PIL	Extensibily index of bread dough	Rad	Radioactivity
	DP Drgestibihty of proteins	PC	Protem content	RH	Relatrve humidity
	EA Ezymatic activity	PIns	Presence of msects	S	Smell of grams
	GH Gram hardness	PMIt	Presence of mites	SW	Specific weight
	GV Gerrrunation vigour	PmO	Proliferation of rmcro-orgarusms	Va	Vanety
	HF HIdden form of msects	PR	Pesticide residues	W	Bakig index
	lIM Heavy metals	Q	Gram quality	W1000	Weight of 1000 grams
			QIpC	Intrmsic & physico-chemical quality	YH	Year of the harvest
				Legend
				o meaningful group o label of the meaningful group
				mes-X measured value of X
				obs-Y ObselVed value of Y CD qualitatIve value of Z
				----..	Influences
				--..	gives-informatIon-on
			o

Table 2 .

 2 The quahtres space QS.

	QuahtIes		Examples of mterpretation
	Vvh	very very high excellent	very cold very very dry
	Vh	very lngh	very good	cold	very dry
	H	lngh	good	fresh	dry
	M	medrum	average	normal	normal
	L	low	mediocre	half-warm	hunud
	VI	very low	bad	warm	very hunud
	Vvl	very very low very bad	very warm	too humid

Table 3 .

 3 Example of projection tables q-sw and q-mc of respectively the specific weight measurement (sw ) and moisture content measurement (me) into the qualities space QS

		SpeCIfIC weight		MOIsture content
		(kg !HI)		(%)
	q-sw	Qualitative value	Weighmg factor	q-mc	Quahtatrve value	Weighmg factor
	sw~80	H	o 9	mc~12	vvh	0.95
	80 > sw ~75	M	0.9	12 < me ~13	vh	0.95
	75 > sw	L	0.9	13 < me ~14	h	0.95
				14 < me ~15	m	0.95
				15 < me ~16	I	0.95
				16 < me	vi	0.95

Table 4 .

 4 Example of combination table (sw-mc) of specific weight (SW) and moisture content (MC) of baking wheat into the qualities space QS; and its associated weighing factors.To describe the grain conditions we used a simple vocabulary, close to the one of store-keepers.

					MC
		sw-mc	vl	L	m	h	vh	vvh
		1	wI	L	I	vI	vI	vI
	SW	m	w I	M	m	m	1	1
		h	wI	H	vh	h	m	I
	wcabulary to describe the grain conditione
							It allows us to
	manipulate concepts such as 'this grain is heated',	'this
	grain is not sound' or 'this grain is too humid' etc. . .. (see
	Table 5).

Table 5 .

 5 Example of vocabulary to describe the grain conditions.

				Me
		sw-mc vI	m
			vvl
	SW	m	vvl In	In
		h	vvl h	vh

ONIC (Dffice National Interprofessionnel des Cles) and ITCF (Jnsutur Technique des Cles et Fourrages) are national orgamsations for gram production and qualrty

Intnnsic & physrco-chemical conditions

Samtary & safety charactensttcs 4 5

Acknowledgements

Part of this project has been financially supported by the Conseil Rional d' Aquitaine.


Thomas Romary
email: thomas.romary@mines-paristech.fr

Fabien Ors

Jacques Rivoirard

Jacques Deraisme

Keywords: Spatial data, clustering, multivariate data, geostatistics

published or not. The documents may come

Unsupervised classification of multivariate geostatistical data: two algorithms Thomas Romary * 1 , Fabien Ors 1 , Jacques Rivoirard 1 , Jacques Deraisme 2

Introduction

In mining assessment, a partitioning of the data is often to be conducted prior to evaluate the reserves. This is necessary to design the mineralization enveloppes corresponding to several lithofacies where the grades of the ore to be mined may have different spatial behavior, in terms of mean, variability and spatial structure. In remote sensing of environmental variables, a similar problem may be encountered when the characteristics of the variable of interest is governed by a hidden variable, e.g. the component of a mixture model, accounting for a particular local behaviour. A typical example in soil sciences consists in the the retrieval of soil classes over a region from the observation of continuous variables.

A natural solution to this problem is to cluster the data. Clustering a dataset consists in partitioning the observations into subsets (called clusters) so that observations in the same cluster are similar in some sense.

Clustering is used in many fields, including machine learning, data mining, pattern recognition, image analysis, information retrieval and bioinformatics [START_REF] Hastie | The elements of statistical learning[END_REF]. It is an unsupervised classification problem where the goal is to determine a structure among the data, with no response variable to lead the process.

While a wide range of methods exist for independent [START_REF] Hastie | The elements of statistical learning[END_REF] or gridded spatial observations (in the image processing litterature), not much attention has been paid to the case of irregularly spaced data. Indeed, in a geostatistical context, one expects to obtain a classification of the data that presents some spatial continuity. This is especially the case with mining data, where the geologist wishes to delineate homogeneous areas in a deposit to facilitate its evaluation and exploitation.

Clustering in a spatial framework has been mainly studied in the image analysis context where the data is organized on a grid. The model is usually a hidden Markov random field. In this model, label properties and pixel values need only to be conditioned on nearest neighbors instead of on all pixels of the map, see e.g. [START_REF] Guyon | Random fields on a network[END_REF] for a review and [START_REF] Celeux | Em procedures using mean fieldlike approximations for markov model-based image segmentation[END_REF] for more recent developments. In [START_REF] Ambroise | Clustering of spatial data by the EM algorithm[END_REF], the authors proposed to use this approach directly to irregularly sampled data using a neighborhood defined by the Delaunay graph of the data. As the length of the edges of the graph are not accounted for in the approach, this neighborhood structure does not reflect a structure in the data, rather a structure in the sampling scheme. This disqualifies this approach especially for mining data, where the samples are located along drillholes: two neighbors on a same drillhole are distant a few centimeters while two neighbors from two different drillholes may be distant several decimeters. [START_REF] Oliver | A geostatistical basis for spatial weighting in multivariate classification[END_REF] were the first to propose a method for the clustering of multivariate non-lattice data. They proposed to modify the dissimilarity matrix of the data, used e.g. in a hierarchical algorithm, by multiplying it by a variogram matrix. This terms to smooth the dissimilarity matrix for close pairs of points. However, this will not enforce the connexity of the resulting clusters, it will rather blur the borders between geologically different areas, making them difficult to differentiate, as our practice showed.

In [START_REF] Allard | Clustering geostatistical data[END_REF], the authors proposed a clustering method based on a mixture of random functions models where an approximation of the expectation-maximization (EM, see [START_REF] Dempster | Maximum likelihood from incomplete data via EM algorithm (with discussion)[END_REF] algorithm is used to estimate the parameters and the labels. It has been later extended to multivariate data in [START_REF] Guillot | Inference of a hidden spatial tessellation from multivariate data: application to the delineation of homogeneous regions in an agricultural field[END_REF]. However this method relies on strong assumptions that are not likely to be encountered in practice: the data are assumed to be Gaussian and data belonging to different clusters are assumed independent. Moreover, the estimation algorithm requires the computation of the maximum likelihood estimator of the random function model at each iteration of the EM algorithm, which involves the inversion of the covariance matrix and is not computationally tractable for large, multivariate datasets. Indeed, a single iteration requires several inversions of a (n×p)×(n×p) matrix, where n is the number of data and p is the number of variables. Using composite likelihood techniques [START_REF] Varin | An overview of composite likelihood methods[END_REF] could be useful to alleviate the computational burden but it will add a degree of approximation while still not allowing to deal with categorical data.

The approaches developped in this paper are model free and do not involve complex computations. Therefore, they are able to process large, multivariate datasets. The first one, already outlined in [START_REF] Romary | Domaining by clustering multivariate geostatistical data[END_REF], is based on an agglomerative hierarchical algorithm with complete linkage (see e.g. [START_REF] Hastie | The elements of statistical learning[END_REF], where the connexity of the resulting clusters is enforced through the use of a graph structuring the data. It only involves the computation of distances along the edges of the graph which has a sparse structure. Its sequential nature makes it practical for reasonable volumes of data. An alternative for large datasets consists however in running first the algorithm on a subsample, then training a supervised classifier and finally applying it to the rest of the data. The second proposed algorithm provides a non-hierarchical alternative to partition the same graph. It is an adap-tation of the spectral clustering algorithm [START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF][START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF] to geostatistical data. The computations involve only sparse matrices, therefore this second algorithm is adapted to large volumes of data.

The paper is organized as follows: in section 2, we describe both algorithms as well as a method to classify newly available data based on the results of a preceding clustering. In section 3, we show the performance of each algorithm on a synthetic dataset as well as on a mining dataset.

Algorithms

Both algorithms proposed rely on the same basic idea. The latter consists in structuring the available data in a graph in the geographical space made of a unique connex component. This graph is then partitioned into clusters either hierarchically or directly by decomposition. The structure thus imposed ensures the spatial coherency of the resulting clusters.

We consider a sample of georeferenced data (x 1 , . . . , x n ) ∈ R n×p , where p is the number of variables, coordinates included. We also consider that the data have been standardized preliminary to the application of the clustering algorithms. It may also be useful to gaussianize the variables, e.g. by anamorphosis [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty, 2nd Edition[END_REF], for skewed data. This preliminary processing allows to make the variables comparable. We describe in this section the different ingredients required to implement both algorithms as well as their core.

Structuring the data

Being either regular or not, the spatial sampling of a geostatistical dataset defines a geometric set, namely a set of points in the geographical space. From this set, a neighborhood system can be built. This can be represented by an undirected graph where each vertex represents an observation and each edge shows the relation of neighborhood shared by close points [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF]. We call this graph the sampling graph. Several methods can be applied to build it such as Delaunay triangulation, Gabriel graph or a graph based on the neighborhood structure defined by moving neighborhood algorithms used in kriging, for instance based on octants (see e.g. [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty, 2nd Edition[END_REF]. Particular shapes can also be obtained by using non-isotropic distances or coordinates transformation. The graph should be parsimonious whilst containing enough edges to support a variety of configurations for the clusters. In our experience, the Delaunay graph and a graph based on a neighborhood selection algorithm give good results. Once the graph G has been built, two observations x i and x j , i = j, are said to be connected if their exists an edge in G linking x i and x j . This is denoted by x i ↔ x j . G can also be represented by its adjacency matrix with general term (G ij ) i,j∈{1,...,n} :

G ij = 1 if x i ↔ x j , 0 otherwise.
Note that an individual is not considered to be connected with itself.

Choice of a distance

The second basic ingredient of both algorithms is a distance or metric measuring the dissimilarity between two observations. The aim of clustering algorithms is to group similar observations, hence the need to define similar. We define the distance d between two observations x i and x j by:

d(x i , x j ) = p k=1 p l=1 ω k,l d (k,l) (x (k) i , x (l) j ), (1)
where (ω k,l ) (k,l)=(1,...,p) 2 are the entries of a positive definite matrix Ω and (d (k,l) ) k=1,...,p is a set of coupled distances, each one chosen according to the corresponding couple of variables. d is therefore a weighted sum of distances. The weights are to be chosen by the user, depending on the relative importance the variables should have and their possible correlation.

As noted above, the variables have been preliminary standardized so as to avoid any scale effect between the variables. In practice, Ω is generally chosen to be diagonal and only individual distances are thus involved. The use of the squared Mahalanobis distance, where Ω is the inverse of the empirical covariance matrix could be considered so as to account for possible correlations between variables, but has not proven useful in our experiments.

The individual distances are chosen according to the associated variable: if the latter is quantitative, the squared euclidean distance is advocated from its strong relation with the variogram as a measure of the local continuity; if it is a categorical variable, an ad-hoc distance is used. Such a distance may take the value 0 when both observations have an equal value for this variable and 1 otherwise. Other options are also available, see e.g. [START_REF] Hastie | The elements of statistical learning[END_REF] for a comprehensive view.

It is worth noting that the coordinates are also included in (1). Indeed, although the spatial location of the data is already accounted for by the graph structure, this allows to account for the length of the edges. By doing this, we promote short connections.

Concerning the setting up of the weights, we generally recommend to put 5% to 30% of the total on the coordinates and to set the other variables to 1 at a first guess, then to progressively tune the weights of the variables according to the outcome of the algorithm.

Geostatistical Hierarchical Clustering

The distance defined above is only valid between pairs of observations.

Agglomerative hierarchical clustering algorithms require a linkage criterion which specifies the dissimilarity of sets as a function of the pairwise distances of observations in the sets. Lance and Williams formula [START_REF] Lance | A generalized sorting strategy for computer classifications[END_REF] enables the use of a unique recurrence formula to update the distances when merging two clusters for a large family of criteria, including the maximum, minimum or average distance, respectively named complete, single or average linkage criteria or Ward's criterion which computes the intra-cluster variance, see e.g. [START_REF] Milligan | Ultrametric hierarchical clustering algorithms[END_REF].

In our context, the spatial continuity needs to be taken into account during the linkage process. In the proposed algorithm, two clusters can merge if and only if they are connected in the graph structure G. When two clusters merge, the resulting cluster inherits all connections of its components. This point is the only departure from the original hierarchical clustering algorithm.

The geostatistical hierarchical clustering algorithm (GHC) is described in pseudo code in algorithm 1 under the complete linkage criterion.

Algorithm 1 Geostatistical Hierarchical Clustering algorithm (GHC)

1: Compute the distance matrix D ∈ R n×n , such that D ij = d(x i , x j ), j < i, if i ↔ j, 0 otherwise 2: repeat 3:
Find k and l, k < l, such that D lk = min {i,j,i↔j} D ij

4:

Merge k and l in {kl}, and update D such that

D ki = max(D ki , D li ) if i ↔ {kl} and i < k D ik = max(D ik , D li ) if i ↔ {kl} and k < i < l D ik = max(D ik , D il ) if i ↔ {kl} and i > l
notation i ↔ {kl} means i is connected with the cluster {kl}, that is i ↔ k or i ↔ l.
Since two clusters are merged when they realize the minimum distance among the connected pairs of clusters, they may not realize the minimum distance in absolute, depending on the chosen linkage criterion. In particular, more dissimilar points may merge into clusters before having merged points which are actually more similar but not directly connected. That is why we advocate the use of the complete linkage criterion which is, to our knowledge, the only way to preserve the ultrametric property in our algorithm. The ultrametric property means a monotonic increase of the dissimilarity value of the clusters, see [START_REF] Milligan | Ultrametric hierarchical clustering algorithms[END_REF] for further details. In particular, the ultrametric property allows to build a dendrogram.

The dendrogram is a very practical tool to select the final number of clusters, see an example in figure 1. It represents the evolution of the intra-cluster dissimilarity along the agglomeration process. A long branch means that the merge between two clusters leads to a much less homogeneous one. Therefore the tree should be pruned at the level where the branches are long. The number of pruned branches gives the number of clusters to consider, 2 in the example of figure 1.

The computational efficiency of this algorithm relies on the graph structure employed and especially on the number of connections. Indeed, only the distances between connected points are required at the beginning of the algorithm, which makes the matrix D sparse and allows fast computations.

Then, the computation of the distances between connected points required at step 4 can be performed on the fly.

Geostatistical Spectral Clustering

The Geostatistical Spectral Clustering (GSC) is an adaptation of the spectral clustering algorithm where the graph used is the sampling graph defined above instead of a graph based on the similarity. Contrarily to GHC, it requires a preselection of the number K of desired clusters and does not rely on an iterative procedure. This is not a major drawback however. Once computed the quantities required for a given maximum number of classes, it is straightforward to compute the outcome for a smaller number of classes.

The different steps of the algorithm are described in algorithm 2.

Algorithm 2 Geostatistical Spectral Clustering algorithm (GSC)

1: Compute the similarity or weighted adjacency matrix W :

W ij = exp - d(x i ,x j ) σ 2 if i ↔ j 0 otherwise (2)
2: Compute the degree matrix D:

D ii = n j=1
W ij 3: Compute the graph Laplacian matrix

L = D -1/2 W D -1/2
4: Compute the K largest eigenvalues of L and form the matrix V ∈ R n×K whose columns are the associated K first eigenvectors of L 5: Apply the K-means algorithm to the lines of V 6: Assign observation x i to the same class the line i of V has been assigned

This algorithm consists in representing the data into an infinite dimensional space (the reproducing kernel Hilbert space associated to the kernel used in (2), here the Gaussian (or radial basis function kernel) where they are easily clustered through K-means.

The parameter σ 2 is chosen as the empirical variance of the variable, following von [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF]. Note that a local adaptive approach could be considered for the setting of σ 2 , as proposed in Zelnik- [START_REF] Zelnik-Manor | Self-tuning spectral clustering[END_REF]. However, this refinement has not proven useful in our practice.

Also, the lines of V can be optionally normalized prior to step 5, as proposed in [START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF]. The differences when using the normalization or not did not appear sensible in our experimentations.

The number of clusters to consider can be chosen by studying the eigenvalues of L. A small eigenvalue signifies that the associated eigenvector is not relevant to discriminates the data. In practice, we advocate to compute a given maximum number of eigenvalues (10 to 20), which corresponds to the maximum number of clusters we want, and then to plot them. A large difference between two eigenvalues means that the smaller one is not so relevant.

As the graph structure is sparse, all the computations required in algorithm 2 can be carried out using sparse linear matrix algebra, which makes GSC computationally efficient and adapted to large multivariate datasets.

Classifying new data

Sometimes, the sampling of the variables of interest on a domain can be performed in several steps. For instance, new drilholes can be added to an initial sampling campaign. In the case where a clustering has already been performed, we may want to classify the new data into the classes resulting from that previous run. An other occurrence when we want to classify data upon the results of a previous clustering is when dealing with very large datasets with the GHC. In that case, we propose to run first the algorithm on a subsample, then train a supervised classifier and finally apply the latter to the remaining data.

It is particularly difficult to incorporate new data into the clustering results

with simple rules. Indeed, when new data are added, the sampling graph gets modified and the outcome of GHC and GSC may change dramatically.

Therefore, the idea developed here is to learn a classification rule based on the initial clustering results. This can be achieved for instance through support vector machines (SVM, see [START_REF] Hastie | The elements of statistical learning[END_REF]). In the case of two classes, the basic principle is to find f (x) = α 0 + N i=1 α i Φ(x, x i ), where the (α i ) i=0,...,N are scalars and Φ a given kernel function, that minimizes as a function of (α i ) i=0,...,N and where the underscript + means the maximum between 0 and the quantity between parenthesis, and λ is a penalty parameter. For multi-class classification, several options are available among which we retain the standard "one versus all" implemented in LIBSVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF]. The penalty parameter λ is set through crossvalidation. Applying the rule to a new observation allows to assign it to an existing class.

N i=1 (1 -y i f (x i )) + + λα t Φα, (3)

Results

Toy dataset

Here, we describe a 2D example on which we have evaluated the performances of several methods including GHC and GSC. We consider a random function on the unit square which is made of a Gaussian random function with mean 2 and a cubic covariance with range 0.3 and sill 1 on the disk of radius 0.3 and center (0.5,0.5) and a Gaussian random function with mean 0 and an exponential covariance with range 0.1 and sill 1 elsewhere. This a class delineating the disk fairly well. It seems that this algorithm tends to generate more compact subsets of the sampling graph.

Each of the five algorithms are applied to 100 realizations of the same random function model, each with a different uniform random sampling. Then we compute the mean, median and 90% percentile of the rate of correctly classified points. Results are summarized in table 1.

GHC exhibits the best performances overall with 85% correctly classified points in average while K-means providing similar results in average, GSC performing the worst with HC and O & W in between. If we look at the median however, GHC has the greatest one with a larger margin. The 90% percentile indicates that in the 10% most favorables cases, GHC misclassified only 0.02% of the points, while all the other algorithms perform worse.

It can also be seen that the 90% percentile are similar for the K-means and the HC. This means that the HC, and GHC (its worse result in this task was a misclassification of almost 50%, seemingly due to a high sensitivity to large values), can sometimes perform really bad, whereas the K-means algorithm gives more stable results, being less sensitive to extreme values.

Indeed, in the presence of very large or very low value, it occurs that the algorithm comes out with a class made of a single point while the other contains all the other observations. In the favorable cases however, HC algorithm works as well as the K-means, while GHC outperforms clearly all other algorithms. Concerning GSC, the results obtained are extremely poor It is worth noting that the drawbacks exhibited by GHC and GSC are far from being prohibitive in practice. Indeed, when applying clustering algorithms to real data the user generally observes the outcome for several numbers of classes. This can be performed easily with both algorithms with a negligible computational cost.

Mining data example

In this section, we present an application of both geostatistical clustering algorithms to an ore deposit. We describe the different steps and exhibit some results.

The first step is to select the data that will be used for the clustering. The following variables are chosen:

• coordinates, X, Y and Z,

• ore grades,

• a geological factor describing the basement vs. a sandy part on top of it,

• the hematization degree.

This choice is made upon an exploratory analysis of the data and discussions with geologists. Some transformations of the data are preliminary performed:

• coordinates are standardized,

• ore grades are log-transformed and standardized,

• the degree of hematization is transformed into a continuous variable, then standardized.

The next step consists of building the sampling graph connections between geographically close samples. The graph is here built from the neighbouring structure induced by the moving neighbourhood kriging algorithm of Isatis R 2013 [START_REF] Geovariances | Isatis technical references[END_REF]. At each point, the space is split into 16 hexadecants: 8 above and 8 below the horizon. One neighbor per hexadecant is authorized at most for each point with no more than 2 from the same drillhole. The search ellipse is of infinite size so as to connect even possibly distant points. The angles of the search ellipse are chosen so that to take into account the horizontal shape of the mineralization of the deposit.

Then the dissimilarity matrix is built. All variables listed above are used.

A particular distance for the geological factor is considered: it is chosen to be 1 when the samples have different factor values and 0 otherwise. This distance is scaled to maintain the coherency with the other individual distances. Weights are set step by step, as advocated in section 2.2: we begin by giving an equal weight to all variables with a 30% contribution to the coordinates. Finally, the contribution of the coordinates is lowered to 10% while the other variables are assigned equal weights. The same set of weights is used for both algorithms. Practice shows indeed that setting low weights to the coordinates leads to better results, as the spatial aspect is already somehow taken into account by the sampling graph. However, the coordinates needs to be included in the distance so as to account for different length of the edges in the graph. This is especially important for drillholes data where two neighbors along a drillhole are generally much closer than two neighbors belonging to two different drillholes.

Finally, we can run both GHC and GSC algorithms described in section 2.

We choose to represent 6 clusters as the intra cluster dissimilarity at that step of the GHC shows a great increase. The results are depicted in figure 5 for GHC and 6 for GSC.

GHC separates the basement into two classes, the black one being richer than the red one. Note that the black cluster is mainly present in the middle of the deposit. The sandy part on top of the basement is splitted into 3 separate classes plus one single observation (in cyan), see figure 5. The discrimination between the three sandy classes seems to rely on geographical considerations.

As for GSC, it splits the basement into 5 classes and puts every observation on top of it into one single class. Some similarities can be observed between the clustering results obtained with the two algorithms however. In particular, both make a clear distinction between the basement and the sand on top of it, emphasizing the variable 'geology'. They also both exhibit the desired connexity properties. Both also reveal a high grades area in the center of the deposit (the black cluster in both figures), whose retrieval was the goal of the experimentation. As already noticed in the previous paragraph, GSC tends to produce more compact clusters than GHC who can follow awkward routes along the graph.

-2 0 1 2 -2 X -2 0 1 2 -3 -1 1 3 -2 0 2 -4 -2 0 -1.5 0.0 1.5 -2 0 1 2 Y -2 0 1 2 Z -3 -1 1 3 U_ppm -2 0 2 geology -4 -2 0 -1.5 0.0 1.5 -1.5 0.0 1.5 degree_hematite a. 1 Y 0 -1 -2 -1 0 1 2 -3 -2 -1 0 1 2 Z X 3 b.

Conclusion

In this paper, we presented two clustering procedures adapted to irreg- are grateful to all members of the consortium for helpful discussions, namely Anglogold Ashanti, Areva, BHP Billiton, Eramet and Vale.

Figure 1 :

 1 Figure 1: Example of a dendrogram

Figure 2 :

 2 Figure 2: One realization of the random function a. and sampling performed b.

  Figures 3 and 4 show the results obtained by each five methods on the realization depicted in figure 2. Each subpicture represents the dataset on scatterplots with respect to the coordinates (x and y) and the sampled value (Z). K-means (a.) identifies well the central area. The result lacks of connexity however. In particular, large values outside of the disk are classified as belonging to the disk and low values within the disk are missclassified as well. It can be seen that the method only discriminates between low and high values of Z: the limiting value between the two clusters can be read as more or less 0.5. HC (b.) also discriminates between low and high value but the limiting value is higher, around 2. To sum up, those two classical methods in an independent observations context fail to produce spatially connected clusters. O & W's approach has been tested with various variograms and variogram parameter values but it never showed any structured result (c.). Our interpretation is that multiplying the dissimilarity matrix by a variogram may erase some dissimilarities, inducing a loss in the structure of the data. The GHC algorithm succeeds in providing a clustering with spatial connexity (d.) though non perfect. A part of the area surrounding the disk is misclassified however. If we turn back to the complete realization in figure 2 a, we can see that the misclassified area corresponds to high values of the realization around the border of the disk that are very close to the values taken inside the disk and are thus difficult to classify correctly. Finally, the GSC algorithm performed a geometrical classification by making a cut along the axis of the first coordinate, see figure 4. However, when looking at the result obtained when asking for five classes, it provided

Figure 3 :

 3 Figure 3: Results of K-means a., hierarchical clustering b., Oliver and Webster's method c. and geostatistical hierarchical clustering d.

Figure 4 :

 4 Figure 4: Results of GSC for two a. and five classes b.

Figure 5 :

 5 Figure 5: Resulting clusters for the GHC algorithm from the variables point of view a. and in 3D b.

Figure 6 :

 6 Figure 6: Resulting clusters for the GSC algorithm from the variables point of view a. and in 3D b.

  ularly sampled spatial data. Both algorithms allow to process large multivariate datasets. They rely on a partition of a graph structuring the data in the geopraphical space, thus ensuring the spatial coherency of the resulting clusters. Two applications have been provided, the first one on a toy example and the second on mining data.The results shown on the toy example validate both algorithms as they are able to produce compact, connected clusters. The results obtained for the mining application are also satisfactory as they highlight a homogeneous area with high grades. Thanks to the sequential nature of GHC, it generates a whole ensemble of coherent clusterings that can be useful to the user: he can visualize the results at different hierarchical levels which helps the interpretation and the choice of the final number of clusters for the end user. Note that GSC does not enjoy this property as the results may change dramatically from one desired number of clusters to another. The main drawback of GHC is its limitation to datasets of reasonable size. It becomes slow when the number of observations goes beyond 10000. In the case of large datasets, a two step approach based on subsampling and supervised classification is proposed.Finally, setting the distance used to compute the graph and the weights associated to each variable allows the practitioner to get different clusterings, according to its knowledge of the geology and the variables he wants to be emphasized in the results. The main difficulty in handling these algorithms is their sensitivity to the different parameters used. Moreover the results are difficult to validate except from the computation of indices of compactness of the clusters or of heterogeneity between them. They are mostly to be validated by the eye of the practitioner whose knowledge of the data should guide in the step by step parameterization of the approach.

Acknowledgement

This work has been partly funded by the Geological and Geostatistical Domaining (G 2 DC) consortium, conducted by Geovariances. The authors