Frédéric Clerc
email: frederic.clerc@inrs.fr

Raymond Vincent

Assessment of occupational exposure to chemicals by air sampling for comparison with limit values: the influence of sampling strategy

This study was aimed at quantifying the impact of sampling duration and the number of measurements taken on the quality of assessing occupational exposure to toluene. To this end, a measurement database was built, based on four campaigns carried out in an industrial printing facility. Five Homogeneous Exposure Groups (HEG) were set up and between 120 and 290 individual measurements lasting from 2 to 8 hours were collected for each of them. These measurements were performed with the objective of comparing them to the 8-hour Occupational Exposure Limit (8h-OEL). The resulting data was used to define a reference exposure profile per HEG: the "gold standard". This exposure profile corresponds to a log normal distribution of measurements from which compliance/non-compliance with the 8h-OEL decision is derived. To simulate the possible sampling strategies used by industrial hygienists, six scenarios were defined, each containing a different number of measurements: two, three, four, six, nine and twelve measurements performed per HEG, over different working days and different seasons of the year. The measurement values per scenario were simulated by sampling from the real measurement per HEG. For each scenario, 1,000 simulated exposure profiles and corresponding simulated compliance decisions were computed. They were compared to the gold standard compliance decision using statistical indicators. Three methods were used for computing the simulated compliance decision: 1) the 95 th percentile must be lower than the 8h-OEL; 2) the exceedance fraction with respect to the 8h-OEL must be lower than 0.1% (as defined by standard CEN 689, Appendix D); 3) the 70% Upper Confidence Limit (70%UCL) of the exceedance fraction with respect to the 8h-OEL must be lower than 5% (as defined by French regulations). The results show that exposure assessment quality increases with both the number of measurements and sampling duration when using the 95 th percentile and exposure assessment based on French regulations, whereas it decreases when using the standard. Moreover, guidelines for the efficient evaluation of chemical exposure in the workplace can be drawn up to help professional occupational hygienists. Indeed, boundaries can be recommended regarding the number of measurements and sampling duration necessary to obtain a reliable exposure assessment while minimizing effort devoted to sampling and analysis.

INTRODUCTION

Assessing the occupational risks linked to the chemical agents present in working atmospheres is based on two main parameters: the intrinsic danger of the agent in question and the exposure of workers to this agent. An agent can take different exposure routes: ingestion, cutaneous contact and inhalation. This paper covers only exposure by inhalation of an agent. To protect exposed employees, the 8-hour Occupational Exposure Limit (8h-OEL) values have been defined on the basis of the toxicity specific to each chemical agent. These values are set out in regulations at both European and French levels. Recently, the transposition of several European directives (EU-OSHA, 2000) (EU-OSHA, 2004) (EU-OSHA, 2006) (EU-OSHA, 2009) has resulted in the publication of strict or indicative regulatory Occupational Exposure Limits (OELs) for over one hundred chemical agents. There are two types of 8h-OEL: short-term exposure limits, which are intended to protect workers from the acute toxicity of a chemical agent, and long-term exposure limits TLV-TWA 8h aimed at protecting workers from the effects of long-term toxicity.

To compare the exposure of workers to an 8h-OEL, a sampling strategy must be used (Leidel et al., 1977) [START_REF] Rappaport | An Exposure Assessment Strategy Accounting for Within-and Between-Worker Sources of Variability[END_REF] (AFNOR, 1995). The principle consists in defining homogeneous exposure groups (HEGs), namely groups of workers carrying out similar tasks and who are considered to be similarly exposed to a given pollutant. Samples are taken from among the workers of the different HEGs formed. The different concentration measurements make it possible to build an exposure profile. In order to assess the exposure profile with respect to the 8h-OEL, various methods can be applied, such as directly comparing the measurements with the 8h-OEL or using simple indicators like averages. The common approach is to assume that a worker should not be exposed to high concentrations for more than 5% of their working time, taking into account any uncertainty on the measurements. To do this, statistical methods are widely used, with the hypothesis of a lognormal distribution of the measurements. The log normal distribution is defined by two parameters: the geometric mean (GM) and the geometric standard deviation (GSD). Based on the log normal distribution, indicators intended to quantify the exposure of the workers of the HEG can be calculated. In particular, the 8h-OEL exceedance fraction -which must not exceed 5% -and the 95 th percentile of the distribution are often used. For each of these statistical indicators, confidence intervals can be derived, and the compliance assessment can be based on the Upper Confidence Limits (UCL) of these intervals [START_REF] Milz | Descriptive Statistics, Inferential Statistics and Goodness-of-Fit. In Ignacio JS and Bullock WH[END_REF] (British Occupational Hygiene Society, 2011).

Assessing the exposure level per HEG is a complex issue as the measurements are subject to considerable variability depending on the emissions generated by the process, the tasks performed by the operator, the season or the type of protective systems employed [START_REF] Rappaport | Assessment of Long-Term Exposures to Toxic Substances in Air[END_REF]. Moreover, sampling implies spatial and temporal variability within an HEG depending on the individual and their exact activity within the HEG [START_REF] Rappaport | An Exposure Assessment Strategy Accounting for Within-and Between-Worker Sources of Variability[END_REF] [START_REF] Symanski | A Meta-Analytic Approach for Characterizing the Within-Worker and Between-Worker Sources of Variability in Occupational Exposure[END_REF]. This variability leads to a high degree of uncertainty in the exposure assessment, resulting in the need to take the most representative air samples of the real working conditions and tasks carried out by the operator throughout their shift. This means taking samples whose duration is appropriate to the work and in sufficient numbers to ensure compliance with 8h-OEL. The problem of optimising sampling strategies has been examined from different viewpoints. In particular, two studies tackled the problem of optimizing sampling strategy by estimating the appropriate balance between the number of measurement days versus the number of workers [START_REF] Lampa | Optimizing Occupational Exposure Measurement Strategies When Estimating the Log-Scale Arithmetic Mean Value-An Example from the Reinforced Plastics Industry[END_REF] [START_REF] Chen | Sampling Strategies for Occupational Exposure Assessment under Generalized Linear Model[END_REF].

This study concerns the relevance of the choice of the duration and numbers of samples required to ensure the reliability of the exposure assessment. Different scenarios (combinations of number of samples / durations) are explored. For each scenario, simulations are performed and an average "compliance" or "non-compliance" decision is compared to a reference "gold standard" decision.

Three methods are used for establishing the "compliance" or "non-compliance" decision. The first, compares the 95th percentile of the probability distribution to the 8h-OEL, while the second uses the CEN 689 -Appendix D standard in which the 99.9 th percentile is compared to the 8h-OEL (i.e. exceedance fraction is compared to a threshold of 0.1%). Finally, the third uses French regulatory texts in which the 70% UCL of the exceedance fraction is compared to a threshold of 5%. The study is based on a case studied in the industrial printing sector where workers are exposed to the toluene contained in the inks.

METHODS

Industrial printing

The exposure measurements were performed in an industrial printing plant using the rotogravure technique. Two workshops were investigated: the first housed three identical production lines, while the second housed only one specialised line set aside for pre-print tests. Each machine uses eight cylinders 3.60 m in width for double-sided, four-colour printing. The machines operate 24 hours a day, with three teams working shifts on the production lines and only one team working on the test line. The average printing time is five hours a shift, depending on the changes made to the print runs or on the problems encountered. Over 150,000,000 kilograms of paper are consumed every year, while between 120,000 and 130,000 kilograms of ink are consumed per week. The ink comprises 50% toluene, which is used as a thinner. The toluene vapours are recovered (95%) by ventilation and a recycling unit comprising a system for trapping toluene vapours on active charcoal. The majority of this toluene is resold, with a small portion being used for cleaning printing line components.

Toluene is a colourless, flammable liquid whose toxicological properties are well known. Acute toxicity is low and chronic toxicity is moderate. Toluene is an R2 category chemical agent classed as toxic to reproduction by the EU (EC regulation no.1272/2008), the 8h-OEL having been lowered to 76.8mg/m 3 and the 15min-OEL to 384 mg/m 3 on 1 January 2012 (Journal Officiel de la République Française, 2012). However, the results of the present study were obtained before the application of this new 8h-OEL. Thus, the 8h-OEL used in this study was the French one in force before January 2012: 192mg/m 3 .

Homogeneous Exposure Groups

The three production lines are strictly identical and none is dedicated to any special process, hence the five HEGs were formed in accordance with this hypothesis. The winder supplies the printing machine with paper. Driver 1 is stationed at the general control panel of the machine, in a ventilated booth of about 20m 2 . Driver 2 is also stationed in this booth. Their role is to check the quality of the print runs and intervene on the printing machine if needed. The receiver is in charge of palletising the printed material. The last HEG concerns two multi-skilled people: the pre-print test drivers who perform all the tasks mentioned above and which depend on the ad hoc requirements generated by the tests to be carried out. There were 9 workers per HEG (3 lines * 3 shifts), but the function of test driver was ensured by 6 workers (2 workers * 3 shifts). The night shift was not sampled. The workers within the HEG are considered strictly equivalent; no information is available on which worker was sampled. Nevertheless, to ensure maximum variability, sampling was performed as much as possible on different workers.

Sampling strategies

Two types of samples were taken simultaneously: a passive sample, the GABIE badge [START_REF] Delcourt | Performance Assessment of a Passive Sampler in Industrial Atmospheres[END_REF], and an active sample: pump and charcoal tube. Throughout their working day, each worker of each HEG was systematically equipped with a badge and a pump, the tube of which was changed every two hours. The time required to change the tubes was not taken into account. Thus, for one shift, three active samples lasting two hours each and one passive sample for a full eight-hour shift were taken. None of them were non-detects; a total of fifteen 2-h active samples and 4 passive samples analyses were missing for non documented reasons. For these situations, the corresponding full-shift measurement was removed. Due to the sequential nature of tube replacement, two types of composite measurements could be performed: two four-hour time periods and one six-hour time period. Thus four types of measurements are studied in the rest of this paper: measurements performed on full-shift passive sampling and the measurements performed on 2-h, 4-h and 6-h active sampling (Figure 1). Note that the two 4h active samples are not strictly independent since the halfway 2h sample is used in both.

Data type

Sampling and analytical technique

The sampling devices (badge or tube) were worn by the workers in the immediate vicinity of the upper airways.

The active sampling device is composed of a tube containing two charcoal stages of 100 mg and 50 mg (SKC®, connected by a flexible tube to a sampling pump ensuring a controlled flow rate of 200 cm 3 /min (Gilian®-LFS 113). The workplace air is drawn through both charcoal stages. The first captures the toluene vapours while the second is used as a control of the saturation of the first. The flow rate of the pump is measured at the start and end of sampling by a soap-bubble electronic flow meter (Gilian-Gilibrator®). Determining the quantity of toluene fixed to the GABIE badge or to the charcoal tubes was done by gas chromatography and internal calibration under the following conditions: (1) Desorption of the toluene vapours using 5 ml of carbon sulphide and (2) analysis of the solution by gas chromatography (CPWAX 52 CB capillary column, L = 50 m, or SPB1 column, L = 60 m, at a temperature of 90 °C) with flame ionisation detection (FID), and internal calibration with nbutyl acetate (NIOSH, 2003).

The GABIE badge takes the form of a plastic casing about 4 cm in diameter containing 550 mg of active charcoal. The surface area of radial diffusion is 7 cm 2 . The calculation of employee exposure to toluene was based on the quantity of toluene captured, the toluene diffusion velocity (36.6 cm 3 /min) and the duration of badge exposure. It was demonstrated that under the same experimental conditions, measurement using the GABIE badge is comparable to that obtained by a tube and pump [START_REF] Delcourt | Performance Assessment of a Passive Sampler in Industrial Atmospheres[END_REF].

Exposure profile and method for compliance decision

The decision to be made based on the series of measurements is either "compliance with 8h-OEL" (the null hypothesis Ho) or "non-compliance with 8h-OEL" (alternative hypothesis Ha). Three methods of obtaining a decision are explored, depending on three descriptors with the assumption of a log-normal distribution of measurements per HEG:

1. The comparison of the 95 th percentile category X 0.95 with the 8h-OEL. The 95 th percentile of the measurements is the concentration value such that 95% of the measurements are below this value (Figure 2, right). Its exact value is classed in five categories depending on the 8h-OEL of toluene, as shown in 2. The comparison of the 99.9 th percentile with the 8h-OEL. This is inspired by the CEN689 Appendix D procedure. If X 0.999 ≤ 8h-OEL then the 8h-OEL has probably been complied with: reject Ha. If X 0.999 > 8h-OEL but X 0.95 ≤ 8h-OEL then no decision can be made: uncertainty. If X 0.95 > 8h-OEL then the 8h-OEL has probably been exceeded: Ha is most likely true. In this study, the uncertain and excess categories are lumped together: if X 0.999 is above the 8h-OEL, then the 8h-OEL is considered exceeded: Ha is most likely true. The exact vocabulary of Appendix D of standard CEN689 (NF, 1995) refers to a threshold of 0.1% and 5% concerning the exceedance fraction. This descriptor can be determined by computation, from the parameters of the lognormal distribution (Figure 2, left) or directly from the series of measurements if the number of values available is large enough.

3.

The French regulatory compliance assessment. As described in the French regulations, sampling strategy and statistical analysis of the exposure measurement results are used for assessing exposure. It is a two stage procedure. Stage 1: a minimum of three measurements per HEG are taken during the initial intervention in an enterprise. If all the measurements are below 10% of the 8h-OEL, then the 8h-OEL is probably complied with: reject Ha. If at least one of the measurements is higher than the 8h-OEL then the 8h-OEL is probably exceeded: Ha is most likely true. Stage 2: If there is no conclusion from stage 1, two other interventions including at least three measurements per HEG are taken in the year at different periods. Every time a measurement exceeds the 8h-OEL, then the 8h-OEL is probably exceeded. When at least 9 measurements have been collected, the 70% upper confidence limit (70%UCL) of the exceedance fraction is calculated. If 70%UCL is lower than 5% (0,05), the 8h-OEL is probably respected: Ha is rejected. Otherwise Ha is most likely true: corrective actions have to be performed (Journal Officiel de la République Française, 2009) (Ogden & Lavoué, 2011).

Simulation of campaign organisation scenarios

For each sampling duration , different scenarios of different numbers of measurements corresponding to different foreseeable situations are formulated. Scenario 1: Two measurements per HEG are made the same day, season indifferent. Scenario 2: Three measurements per HEG are made on three consecutive days, season indifferent. Scenario 3: Four measurements per HEG are made, two measurements per day on two consecutive days, season indifferent. Scenario 4: Six measurements per HEG are made, one measurement per day on three consecutive days over two consecutive seasons. Scenarios 5 and 6: Nine and 12 measurements per HEG are made, one measurement per day on three consecutive days over three and four consecutive seasons. This approach corresponds to the regulatory approach in France.

These scenarios represent the different strategies which could be used by industrial hygienists, depending on the resources available for the exposure assessment.

Once 1,000 runs per scenario have been performed, the corresponding exposure profiles are calculated, averaged and compared to the gold standard exposure profiles. To make this comparison, three metrics are computed, based on the descriptors. These metrics represent the average quality of the "compliant" / "non compliant" decision per scenario, depending on the numerical procedure used:

 The error when using the 95 th percentile category comparison. For each run, the category of the 95 th percentile is computed and a decision is made: category 4 (X 0.95 > OEL) leads to noncompliance, other categories lead to compliance. Each decision is compared to the gold standard decision. For example, if the gold standard category is 2 (X 0.95 < 50%OEL), then a simulation leading to a category 1 (X 0.95 < 10%OEL) is correct because it is the same decision: compliance. The percentage of runs resulting in correct and incorrect decisions is then computed for each scenario.

 The error when using the 99.9 th percentile (CEN689 Appendix D standard). For each run, the decision with regards to the X 0.999 is computed and compared to the gold standard decision. Following the same principle, the percentage of runs resulting in correct and incorrect decisions is then computed for each scenario.

 The error when using the French regulatory compliance assessment. For each run, the simulated decision of compliance or non-compliance is computed and compared to the gold standard decision. When the simulation methodology does not exactly correspond to the specifications (scenarios 1 to 4), the decision from stage 2 is established on the basis of the number of measurements taken, while keeping the conditions of immediate decision from stage 1 active. Again following the same principle, the percentage of runs resulting in correct and incorrect decisions is then computed for each scenario.

RESULTS

Gold standard exposure profile

In Table 2, the exposure profiles and decisions based on all available measurements for each HEG are calculated using the AltrexChimie software tool (INRS, 2010). Table2. Exposure profile parameters. n is the number of samples; "X 0.95 " is the 95 th percentile computed with a log normal hypothesis; 95%CI on the 95 th percentile is the 95% confidence interval on the 95 th percentile; X 0.95 category is the AIHA-defined category; X 0.95 decision is C for compliance and NC for non-compliance; "exceedance fraction" is the exceedance fraction calculated on the basis of a log-normal hypothesis; X 0.999 (CEN689 decision) is C for compliance, U for uncertainty and NC for non-compliance; the 70% Upper Confidence Limit (UCL) on the exceedance fraction is used for the French regulatory compliance decision, it must not be above 5%; C is for compliance, NC is for noncompliance and NC* is non-compliance because at least one measurement exceeds 8h-OEL; "% excess measurements" is the percentage of measurements exceeding the 8h-OEL.

For the same HEG, the decision of compliance or non-compliance differs depending on the sampling duration. Indeed, the measurement variability (measured through GSD) always decreases with the duration of samples whereas GM always increases with the duration of samples. According to a 1% threshold for the p-value of the Shapiro-Wilk's test, the hypothesis of a log-normal distribution is rejected in all cases for the "driver 2" HEG (this is due to the particular fact that 11% of measurements exceed the 8h-OEL). The log normal distributions of 8-h measurements per HEG are presented in Figure 3. It is assumed that the series of full-shift measurements collected is fully representative of the usual measurable toluene concentration in the workplace: around 50 measurements per HEG were performed and variability was very low in all cases. Consequently, the gold standard decision is based on the exact percentage of measurements exceeding the 8h-OEL, which must not exceed the threshold of 5%. The conclusion is that:

 The exposure of the workers of the "winder", "driver 1", "receiver" and "test driver" HEGs is probably lower than the 8h-OEL (reject Ha);

 The exposure of the workers of the "driver 2" HEG is probably higher than the 8h-OEL for toluene (cannot reject Ha).

Figure 1 .

 1 Figure 1. Organisation of measurement strategies: single and composite The average obtained from three 2-hour measurements represents the mean value obtained in the situation where the IH monitored exposure for 2-h, whereas the 6-h measurement represents an almost full-shift. In the remainder of this paper, the series of measurements are identified as: Full-shift -Shift measurement with passive sampler; 2-h -Average value of three consecutive 2-h measurements with active samplers (start, halfway, end of shift); 4-h -Average value of two consecutive composite 2-h active samples (start to half-way of shift and half-way to end of shift) 6-h -Shift measurement obtained with the three consecutive 2-h active samples The overall campaign consisted of four in-situ interventions lasting 4 days each, organised over different periods over the years 2000 and 2001 (spring 2000 -57 full-shift samples, autumn 2000 -55 full-shift samples, winter 2001 -64 full-shift samples, summer 2001-full-shift 69 samples) in order to integrate seasonal variability.

Figure 2 .

 2 Figure 2.Log-normal distributions. On the left, graphical representation of the exceedance fraction. On the right, graphical representation of the 95 th percentile X 0.95

Figure 3 .

 3 Figure 3.Log-normal distributions of all HEGs: gold standard distributions obtained with the full-shift measurement series used to establish the gold standard exposure profiles.

Table 1

 1

	(Mulhausen et

The simulations were performed using algorithms developed with Java 1.7. The probabilistic library used for log-normal computations was JScience 5.0.

When the compliance decision is based on the 95 th percentile category, the higher the number of measurements, the better the accuracy. Whatever the duration, the proportion of correct estimations gradually increases to converge towards a plateau from 8-10 measurements onwards. When 10 or more measurements are performed, the full-shift measurements yield the best results, giving a classification error rate of less than 20%. The 4-h and 6-h measurements differed little, resulting in an error rate of less than 30%. The 2-h measurements clearly resulted in a poor error rate, above 40% (Table 3). Table3. Average proportion of classification errors of the X 0.95 with respect to sampling duration and number of measurements. The columns for 5, 7, 8, 10 and 11 measurements are extrapolated linearly from adjacent values. Values above 50% indicate that the proportion of X 0.95 estimations is more likely to be in the wrong category.

Full

In the simulation with 9 full-shift measurements (20% error rate), the false negative and false positive decision rates per HEG are reported in Table4. A so-called false positive is reported when the correct decision is "compliance" but the computed decision is "non-compliance". A so-called false negative is reported when the correct decision is "non-compliance" but the computed decision is "compliance". The false negative is therefore more problematic [START_REF] Hewett | A 21 st century exposure assessment strategy -Minimizing sampling strategy false negative and false positive decisions[END_REF]. When the compliance decision is based on the threshold of 0.1% of the exceedance fraction; the mean percentage of errors increases as the number of measurements increases. Taking 9 full-shift measurements leads to an error rate of 60%. This also applies to nine 6-h and 4-h measurements, which result in rates of about 60%, while the order of magnitude is 70% for 2-h measurements (Table 5).

HEG

Full-shift 41% 50% 42% 49% 56% 57% 59% 60% 61% 61% 62% In the simulation with 9 full-shift measurements (60% error rate), the false negative and false positive decision rates per HEG are reported in Table 6. Table6. Average false positive and false negative rates when the decision is made on the CEN689 per HEG for scenario 5 (1000 simulated datasets of 9 measurements of 8-h).

HEG

When the quality of exposure assessment is based on French regulations (Journal Officiel de la République Française); the percentage of decision errors decreases as the number of measurements increases. In almost every case, the quality of estimation also benefits from the increase in sampling duration. This is the decision mode that most limits the quantity of errors made compared to the gold standard. Indeed, taking nine full-shift or 6-h samples results in a decision error rate of about 15%, while the order of magnitude is 20% for 4-h samples and 40% for 2-h samples ( In the simulation with 9 full-shift measurements (17% error rate), the false negative and false positive decision rates per HEG are reported in

Discussion

The six scenarios explored in this paper were defined empirically. However, the error values from scenario number 2 diverge slightly from the global trend: the number of errors made is higher. This scenario uses a strategy based on 3 measurements per day over 3 consecutive days; whereas other scenarios always include at least two measurements made the same day on two different workers.

Although this was not verified, the hypothesis is that the amount of variability accounted for by scenario no.2 is higher than in other scenarios, leading to higher GSD and more errors. In other words, scenario no.2 minimizes autocorrelation [START_REF] Symanski | An Investigation of the Dependance of Exposure Variability on the Interval between Measurements[END_REF]) [START_REF] Wackernagel | Applying geostatistics to exposure monitoring data in industrial hygiene[END_REF]. In [START_REF] Lampa | Optimizing Occupational Exposure Measurement Strategies When Estimating the Log-Scale Arithmetic Mean Value-An Example from the Reinforced Plastics Industry[END_REF]) [START_REF] Chen | Sampling Strategies for Occupational Exposure Assessment under Generalized Linear Model[END_REF], it is shown that sampling a large number of workers is preferred to sampling a large number of days, which is different from our observation: in the study presented in this paper, sampling a large number of days seem to increase the error rate.

The "uncertainty zone" is defined in standard CEN689, Appendix D, as the exceedance fraction between 0.1% and 5% (equivalent to percentiles X 0.999 and X 0.95 ). For this study, the choice was made to strictly comply with the terms of the standard, thereby considering that since the "uncertainty zone" is different from "compliance", it must be "non-compliance". The other option consisting in considering the uncertainty zone as "compliance" was not studied but its definition is analogous with the 95 th percentile point estimate. The results observed when computing the error when using X 0.999 (standard CEN689, Appendix D), highlight a large number of errors because most of the simulations led to the uncertainty zone. Concerning the results of table 4 in particular, the exceedance fraction for 3 out of 5 HEGs places them in the "uncertainty zone" (driver1, test driver and receiver). In this case, the "true" decision is "uncertainty", no matter how many measurements are collected. Since "uncertainty" is considered as "non-compliance" in this study, the percentage of errors increases with the number of measurements. This highlights the major weakness of CEN-689.

Considering the impact of sampling duration on the assessment of French regulatory compliance (Table7), there is a clear difference between a 2-h sampling duration and a longer duration. The simulations based on the measurements stemming from the 6-h samples performed better, on average, than those stemming from full-shift samples (whereas the gold standard decision is based on full-shift samples). This is due to the immediate decision mode from stage 1 that uses the fact that a single measurement exceeding the 8h-OEL results in a non-compliance decision. Indeed, only 3 values stemming from 6-h samples actually exceeded the 8h-OEL, while 6 values stemming from fullshift samples exceeded the 8h-OEL.

The dataset used in this study was compiled under real conditions, but it is in fact almost a textbook case since the workers' occupation is strictly defined and the GSDs are very low. The study results do not reflect cases where the occupations of workers from exposure groups are more heterogeneous. Under such conditions, where variability is assumed to be higher, the uncertainty of the compliance decision would also be higher, thus the conclusions would be affected. Moreover, the data available for this study did not allow simulating the situation where 3 measurements are below 10% of the 8h-OEL.

Nevertheless, for a rigorous, full-shift, 8h-OEL compliance assessment, the CEN689-Appendix D procedure is too restrictive due to the use of a very low 8h-OEL compliance threshold: 0.1%. The "95 th percentile" and "French regulatory compliance assessment" decision-making procedures allow making similar conclusions.

 Three or fewer samples are insufficiently reliable to yield a correct judgement.  Samples of 2-h or less are insufficiently reliable to yield a correct judgement.  It does not appear indispensable to take over 10 measurements as the gain in the quality of judgement obtained is low.

Thus, taking between 4 and 9 6-h or full-shift samples provides a relatively good judgement of exposure. French regulations recommend assessing occupational exposure by following a sampling strategy leading to taking 9 full-shift measurements. It has been noted that in this case study, the decision can turn out to be false in 17% of cases. In the case of the HEG defined as non compliant, the rate of false negative decisions reaches 22.5%, while the average rate of false positive decisions for HEGs defined as compliant ranges from 0% to 22.9%. Nevertheless, by following the same sampling strategy, the use of Standard CEN689 results in 60% errors or uncertainties. The use of a criterion based on the 95 th percentile will be erroneous in 20% of cases. Finally, the strategy recommended by French regulations, while not exempt of faults, minimised the number of decision errors in this case study.

Conclusion

This paper described the impact of sampling duration and number of measurements on the quality of estimating worker exposure. The collection of a large number of workplace air samples in an industrial printing plant allowed formulating a "gold standard" which was used as an absolute reference of occupational exposure for five HEGs with low variability. Four HEGs were considered compliant with regards to an 8h-OEL, but the fifth was not. Simulations were carried out to assess average exposure when different sampling number and duration strategies were applied.

With regards to the simulations performed and sampling time, the results show that the quality of exposure estimation is similar for long-term samples (more than 6 hours), regardless of the number of measurements and the type of calculation used for the exposure decision. In contrast, estimating exposure from 2-h samples for comparison with the 8h-OEL is unreliable, and this sampling duration should be abandoned for this specific purpose. In the context of French regulations (9 measurements), 17% of decision errors were made on the basis of long-duration measurements, while 40% of errors were made on the basis of 2-h measurements.

Concerning the number of samples to be taken, the results show that three measurements are insufficient to judge full-shift occupational exposure in a context of low variability. From four measurements onwards, the results gradually improved until reaching a plateau, the inflection taking place between eight and 10 measurements. Even when the resources invested were excellent, with 12 full-shift measurements, the decision error rate in terms of the French regulation is lower but cannot be zero (13% minimum).

For exposure assessment based on a full-shift, when the HEGs are clearly defined with low variability as presented in this study, the strategy set out in the French texts constitutes a compromise providing reasonably reliable exposure judgement while keeping a reasonable number of samples. For a rougher but also faster and less costly exposure assessment under the same conditions, it appears that taking four long-duration measurements spread over two consecutive days is the minimum necessary in order to obtain non-prohibitive error rates (less than 40%). In this case, an indepth examination of the different indicators, like working conditions and measurement series could provide additional information making it possible to judge whether any over-or underestimation of exposure has occurred. In this context, the judgement of the industrial hygienist is available to compensate for any sampling errors brought about by the sampling strategy. Techniques that combine expert judgement and existing measurements have been developed and their use should increase [START_REF] Hewett | Rating Exposure Control Using Bayesian Decision Analysis[END_REF] [START_REF] Sottas | An empirical hierarchical Bayesian unification of occupational exposure assessment methods[END_REF].


Bag-of-Bags of Words model over irregular graph partitions for image retrieval Yi Ren, Aurélie Bugeau, Jenny Benois-Pineau University of Bordeaux, LaBRI, UMR 5800, F-33400 Talence, France Email: {yi.ren, aurelie.bugeau, jenny.benois-pineau}@labri.fr

Abstract-The paper presents a novel approach, named bagof-bags of words (BBoW), to address the problem of Content-Based Image Retrieval (CBIR) from image databases. The proposed bag-of-bags of words model extends the classical bag-ofwords (BoW) model. An image is represented as a graph of local features on a regular grid. Then irregular partitions of images are built using different graph cutting methods. Each graph is then represented by its own signature. Compared to existing methods for image retrieval, such as Spatial Pyramid Matching (SPM), the BBoW model does not assume that similar parts of a scene always appear at the same location in images of the same category. The extension of the proposed model to pyramid gives rise to a method we name irregular pyramid matching. The experiments demonstrate the strength of our method for image retrieval when the partitions are stable across an image category. The experimental results for Caltech101 benchmark show that our method achieves comparative results as SPM, and is globally more stable.

I. INTRODUCTION

Recent methods in Content-Based Image Retrieval (CBIR) mostly rely on the bag-of-visual-words (BoW) model [START_REF] Sivic | Video google: a text retrieval approach to object matching in videos[END_REF]. The idea, borrowed from document processing, is to build a visual codebook from all the feature points in a training image dataset. Each image is then represented by a signature, which is a histogram of quantized visual features-words from the codebook. Image features are thus considered as independent and orderless. The traditional BoW model does not embed spatial layout of local features in image signature. However, those information have shown to be very useful in tasks like image retrieval, image classification, and video indexing. Ren et al. [START_REF] Ren | Learning a classification model for segmentation[END_REF] put forward a concept of grouping pixels into "superpixels". Leibe et al. proposed to adopt codebooks to vote for object position [START_REF] Leibe | Combined object categorization and segmentation with an implicit shape model[END_REF]. Lazebnik et al. [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] partitioned an image into increasingly fine grids and computed histograms for each grid cell. The resulting spatial pyramid matching method (SPM) clearly improves the BoW representation. Nevertheless, this method relies on the assumption that a similar part of a scene generally appears at the same position across different images, which is not always true.

Graphs are versatile tools to conveniently represent patterns in computer vision applications and have been vastly investigated. By representing images with graphs, measuring the similarities between images becomes equivalent to finding similar patterns inside series of attributed (sub)graphs representing them. Duchenne et al. [START_REF] Duchenne | A graph-matching kernel for object categorization[END_REF] introduced an approximate algorithm based on graph-matching kernel for category-level image classification. Gibert et al. [START_REF] Gibert | Graph embedding in vector spaces by node attribute statistics[END_REF] proposed to apply graph embedding in vector spaces by node attribute statistics for classification. Bunke et al. [START_REF] Bunke | Towards the unification of structural and statistical pattern recognition[END_REF] provided an overview of the structural and statistical pattern recognition, and elaborated some of these attempts, such as graph clustering, graph kernels and embedding etc., towards the unification of these two approaches.

In this paper, we present a new approach for Content-Based Image Retrieval (CBIR) that extends the bag-of-words (BoW) model. We aim at embedding color homogeneity and limited spatial information through irregular partitioning of an image into a set of pre-defined number of graphs. Each partition results from the use of graph partitioning methods, the node of the initial graph being positioned on a dense regular grid of pixels. We compare two graph partitioning methods, namely Graph Cuts [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] and Normalized Cuts [START_REF] Shi | Normalized cuts and image segmentation[END_REF]. The BoW approach is next applied to each of graph partition independently. An image is then represented by a set of graph signatures (BoWs), leading to our new representation called bag-of-bags of words (BBoW). As in the spatial pyramid matching approach [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], we also consider a pyramidal representation of our images. BBoW models are therefore computed at different resolutions, i.e. with different number of graphs. More precisely, a different number of (sub)graphs are considered at each level of the pyramid. The comparison of images in a CBIR paradigm is then achieved via comparison of the irregular pyramidal partitions. We call this pyramidal approach Irregular Pyramid Matching (IPM). Hence in this paper, we try to address a challenging question: will an irregular, segmentation-like partition of images outperform a regular partition (SPM)? Intuitively, it is invariant to the rotation and reasonable shift transformations of image plane. Nevertheless, what can be its resistance to the noise and occlusions? How will it compete with SPM when embedded into pyramidal paradigm?

The remainder of the paper is organized as follows. In section II we briefly introduce the notations and prerequisites. The proposed bag-of-bags of words model is discussed in section III and the irregular pyramid matching in section IV. Section V presents the experimental results and we conclude in section VI.

II. TERMINOLOGY

The input database Ω = (I, V ) is composed of N RGB images I = {I 1 , . . . , I N } and of a visual codebook V = {V 1 , . . . , V B } of length B. The codebook is built using kmeans over the SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] descriptor of all the dense grid points from the database. Let us denote by G j = (V j , E j , W j ) an undirected weighted graph constructed on the image I j . The set V of vertices contains a regularly sampled subset of pixels P of the image and at the limit can contain all of them. The graph edges E connect these vertices with a 8connected neighbourhood system. The adjacency matrix W of size |V| × |V| is defined as:

Wpq = wpq if p, q ∈ E 0 if p, q ∈ E,
where w pq represents the edge-based similarity between two vertices p and q. For each image I j in the database, we aim to partition the graph G j into K disjoint subgraphs {g 1 j , . . . , g K j }, such that ∀k = l, g k j ∩ g l j = ∅ and Gj = K k=1 g k j . We denote this K-way partitioning by Γ K j = {g 1 j , . . . , g K j }.

The graph partitioning problem can be recast as a labeling problem. Given a set of vertices V and a finite set of labels L = {1, 2, . . . , K}, for all node p ∈ V, we are looking for the optimal label l p ∈ L, such that the joint labeling L = {l 1 , . . . , l |V| } ∈ L |V| satisfies a specified objective function.

III. BAG-OF-BAGS OF WORDS MODEL

Our BBoW model has been inspired by bag-of-words model [START_REF] Sivic | Video google: a text retrieval approach to object matching in videos[END_REF] for image description and the idea that an irregular partition into a fixed number of graphs is a step in-between an arbitrary regular partition as in [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], and a segmentation of an image plane which can be strongly redundant. The imperfections of this intermediate partition should be compensated by a statistical nature of the original BoW model for each graph in the partition. Furthermore, we embed BBoW into multiresolution schemes as spatial pyramid matching [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] for using image partitions in coarse to fine manner. For each image, construction of a bag-of-bags of words (see figure 1) can be decomposed into four main steps:

1) Select a reduced number of pixels V from the whole image.

A. Graph nodes selection

We adopt a dense sampling strategy, i.e. sample the image at regular grid points. All the areas contribute equally to the final image representation. We choose 8-pixels spacing and 16-pixels patch size for the SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] features, called Bag-of-Features Grid SIFT (BF-GSIFT) [START_REF] Tuytelaars | Dense interest points[END_REF].

B. Construct initial weighted graph

The set of nodes V of our initial graph G contains all the points resulting from the dense sampling. The edges E are obtained by linking these points with a 8-connected neighbourhood system. As mentioned in the introduction, we have tested two different methods to achieve graph partitioning, namely Graph Cuts and Normalized Cuts. For both, the same weights on the edges are used:

wpq = exp -λ Cp -Cq Σ -1 Cp -Cq T , (1)
where Cp accounts for the mean color vector over a n×n patch centred on point p. We use the YUV color space that has independent color channels and allows to better deal with changes in lightning conditions. We denote Y p = Y (x, y) (respectively U p = U (x, y) and V p = V (x, y)) as the color channel value at point p = (x, y) in the whole image coordinate system. The mean color vector becomes Cp = ( Ȳp, Ūp, Vp) where:

Ȳp = 1 n 2 n u=-n n v=-n Y (x + u, y + v).
The covariance matrix

Σ =   σ 2 Y 0 0 0 σ 2 U 0 0 0 σ 2 V  
is diagonal because of channels' independence. As in [START_REF] Rother | Grabcut: Interactive foreground extraction using iterated graph cuts[END_REF], the covariance is defined for each channel as:

σ 2 = (Yp -Yq) 2
where . denotes the expectation over all the edges (p, q) ∈ E of the graph G.

C. Graph Partitioning

Let us now present how the initial graph G is partitioned into K subgraphs. We investigated two different algorithms for that purpose. In the following, the details for both of them: Normalized Cuts and Graph Cuts are given.

1) Partitioning with Normalized Cuts: Normalized Cuts [START_REF] Shi | Normalized cuts and image segmentation[END_REF] is a well known method in data clustering and image analysis. The principle is to divide the graphs into subgraphs such that the similarity between nodes within a subgraph is greater than the similarity between nodes in separated subgraphs. The algorithm then optimizes two criteria measuring the total dissimilarity between the different partitions as well as the similarity within each partition. The Normalized Cuts method can directly be applied to partition a weighted graph G. In practice, we used the code from [START_REF] Shi | Normalized cuts and image segmentation[END_REF], [START_REF] Cour | Normalized cut segmentation code[END_REF].

2) Partitioning with Graph Cuts: Graph Cuts is another very famous method for graph partitioning in the computer vision community. Graph Cuts [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF], [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF] are based on mincut/max-flow optimization and are used to minimize energy functions containing a data term and a smoothness term. In our case, the energy reads:

E(L) = E data (L) + E smooth (L) = p∈V E d (p, l p ) + (p,q)∈E E s (p, q) = p∈S E d (l p ) + (p,q)∈E E s (p, q).
The first term is the data term. It is only applied to seed points p ∈ S, where S = {s 1 , . . . s K } is the set of seeds. Seed points are nodes that we select in the following way. The image is divided into K regular cells, K being the number of partitions we are looking for. The seed points are then the nodes that are Fig. 2: Seed points. In red: the nodes of the graphs. In other colors: the seed points the closest to the barycentres of these cells (figure 2). The goal of the data term is to fix the label of the seed points. The seed points are then used to impose hard constrains on Graph Cuts. For instance, the seed point s i must end up with label i after the energy minimization process. This idea can be formulated as in [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF]:

E d (lp) = 0 if p = si & lp = i ∞ otherwise.
(

) 2
The smoothness term on the other hand, is directly linked to the weights of the graph, as defined in equation ( 1) :

Es(p, q) = wpq pq 2 1δ(lp, lq) .

(

) 3
where δ(.) is the Kronecker function.

D. Signature of subgraphs

As in the standard BoW approach, the signature of each subgraph g j , j = 1 . . . K is a histogram H j . It is obtained by assigning the SIFT descriptor of each of its nodes p ∈ V j to the closest word in the codebook computed on the whole database. As mentioned before, the codebook is computed from all the dense grid points in a training image dataset. In other words, we do not rely on our graph partitions to construct the visual dictionary.

IV. IPMK: IRREGULAR PYRAMID MATCHING KERNEL

In this section, we explain how our bag-of-bags of words model can be used in an image retrieval system to match images.

Our strategy directly follows the one from Spatial Pyramid Matching [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF]. Let us assume that we compute a sequence of subgraphs at different resolutions r = 1 . . . R of an image. For each image I i , we end up with a set H r i = {H r i,1 , . . . , H r i,Kr } of K r histograms at each resolution r. Example of subgraphs and their histograms are visible in figure 3. As with SPM, at each resolution the number of partitions is given by Kr = 2 2r . In order to compare two images I i and I j , we use the histogram intersection function [START_REF] Jou | Efficient matching of large-size histograms[END_REF]:

I H r i,k , H r j,k = B b=1 min H r i,k (b), H r j,k (b) , (4)
where B in the codebook size.

In our method we cannot directly apply this function to match histograms of a pair of images. At the end of the partitions steps, for two images, the two histograms H r i,k and H r j,k may not correspond one to the other. Indeed, after applying a graph partitioning method, the subgraph attributed Fig. 3: A schematic illustration of the bag-of-bags of words representation at each resolution of the pyramid. At level 0, the decomposition consists of just a single label, and the representation is equivalent to the classical bag-of-words. At level 1, the image is subdivided into four partitions, leading four features histograms, and so on.

to label k can be at any position in the image or represent any object of the image. The spatial arrangement is lost if an image in a database undergoes rotation for instance. To overcome this issue, we first need to reorganize the labels of our histograms. For this purpose, we reorganize the labels between the two sets of histograms {H r i,k } k=1,...,Kr and {H r j,k } k=1,...,Kr . We call this step bipartite subgraphs matching.

To rearrange the labels, we are facing an assignment problem which requires the use of combinatorial algorithm. In our method, we rely on the Hungarian algorithm to minimize the discrete optimal transport between the two set of histograms. The cost matrix D r i,j between the pair of images I i and I j at level r reads:

D r i,j =    d 11 i,j . . . d 1Kr i,j
. . .

d Kr 1 i,j . . . d Kr Kr i,j    . (5)
where

d k,l i,j = B b=1 |H r i,k (b) -H r j,l (b)|
is the L 1 distance between H r i,k and H r j,l , which is equivalent to the histogram intersection kernel. The Hungarian algorithm finds the minimum cost assignment by associating each label k in image I i to one label k ′ = f i (k, r) in image I j thus the pairs of histograms to compare are identified. Now that bipartite subgraphs matching has been achieved, we can directly apply the Spatial Pyramid Matching scheme by using the pyramid match kernel:

κ(Ii, Ij) = I H R i,k , H R j,k ′ + R-1 r=0 1 2 R-r I r -I r+1 = 1 2 R I H 0 i,1 , H 0 j,1 + R r=1 1 2 R-r+1 I H r i,k , H r j,k ′ .
The weight 1 2 R-r before the histogram intersection function, which is proportional to number of subgraphs, allows to give more influence to penalize higher resolutions, reflecting the fact that higher levels (i.e. smaller resolutions) localize the features more precisely.

V. EXPERIMENTAL EVALUATION

In this section, we present the experimental evaluation on Caltech101 benchmark [START_REF] Fei-Fei | Learning generative visual models from few training examples: an incremental bayesian approach tested on 101 object categories[END_REF] and the people-playing-musicalinstrument (PPMI) 1 dataset [START_REF] Yao | Grouplet: A structured image representation for recognizing human and object interactions[END_REF].

The Caltech101 dataset includes 101 distinct object such as faces, watches, ants, pianos, etc... and a background category for a total of 102 categories. This benchmark acts as a common standard for comparison of different algorithms without bias due to different datasets. We follow the same setting as in [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF]. We use exactly the same grid points as generated from the code provided by the authors of [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] to build the initial graphs. The dictionary is built over 30 images per class and the retrieval is run on the rest of the images.

The PPMI dataset emphasizes on understanding subtle interactions between humans and objects. For twelve musical instruments, it contains images of people either playing or just holding the instrument. The twelve binary classes are used in our experiments. For each class, 100 normalized PPMI+ images are randomly selected for training and the remaining 100 images for testing.

A. Parameters evaluation

We start by testing the parameters' influence of our method for image retrieval. Two parameters have been evaluated: the patch size n defined in section III-B and the parameter λ from the edges weight (equation 1). The mean average precisions on the two datasets are presented in tables I and II. For efficiency, these experiments were conducted on smaller databases containing only 4 test images per category. These results highlight that these parameters do not have much influence on the quality of results using either Graph Cuts (GC) or Normalized Cuts (NC). For the following experiments, we therefore decided to use n = 5 and λ = 5.

B. Image retrieval evaluation

The performance of our algorithm is measured by computing the mean Average Precision (mAP) for each category, as described in TREC-style evaluation2 . We also add the computation of the mAP and the standard deviation for the whole dataset. We compare with SPM [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] and with the bagof-words approach which simply corresponds to the results at level 0. Tables III and IV present this global mAP for the two datasets at each level of the pyramid. SPM globally outperforms our method, but we outperform BoW (level 0) both with Graph Cuts (GC) and Normalized Cut (NC). When looking closer to the results, we can nevertheless observe that the pyramidal approach is of no need for SPM on the dataset used. Indeed, the results at level 2 are better than the ones with the pyramid approach. These results have been obtained with the codes from the authors of SPM. On the contrary, the pyramidal approach is a good strategy in our case, as it improves the retrieval performance. We now propose to look into more details to the results for certain categories. A comparison of the mAP on sixteen categories from Caltech-101 dataset is presented on figure 4. The result highlights the fact that for a few number of categories SPM has a much higher retrieval performance than our method. These categories explain the importance difference in the global performance on all the database (as described previously). Nevertheless for most of the categories, we are as good or better than SPM with either Graph Cuts or Normalized Cuts. In order to give a deeper understanding on why our method suits for some catagories and not for the others, let us now look at the graph partitions obtained by the graph cutting methods. Figure 5 presents the results of the partitioning at level 2 (16 partitions) with either Normalized Cuts (NC) or Graph Cuts (GC). The first two lines show images from two different categories for which our method clearly outperforms SPM. For the minaret category, IPMK with Normalized Cuts gives the best results. As can be seen on the first row of figure 5, our method is able to represent the minaret with only one or two partitions. The other partitions in the background are stable in the sense that almost the same partitions can be found in each image. The second row presents the Graph Cuts results on four images from the beaver category. Once again, it can be observed that the partitions are stable. The four last rows on the other hand represents two categories for which our

VI. CONCLUSION AND PERSPECTIVES

In this paper, the bag-of-bags of words model over irregular image partition with graph cut methods and its pyramidal extension, the irregular pyramid matching kernel have been introduced. Graph cuts methods are used to separate an image graph which is a dense grid of points into several subgraphs, each then being described by a bag-of-words model. Our approach incorporates color and limited spatial information into image representation for image retrieval. The obtained results are encouraging, especially when the partitions obtained are stable across all the images from the same category. Nevertheless, this stability is not always insured with this preliminary work. We will now be focusing on improving the weights of the graph in order to reach more stability.

2 )

 2 Build an initial graph G. 3) Partition the graph G into K subgraphs. 4) Compute a signature for each subgraph.The signature of a subgraph is a histogram obtained by assigning each feature of this subgraph to the closest visual word in the codebook. Hence, an image being composed of K subgraphs is characterized by a set of K histograms. In the following, we detail each step of the method.

Fig. 1 :

 1 Fig. 1: Bag-of-bags of words pipeline

Fig. 4 :

 4 Fig. 4: The mean Average Precision for 16 typical categories in Caltech101 method fails. The partitions obtained almost do not exhibit any consistency across the images from the same category. The resulting histograms are therefore each time very different and we can obviously not expect any good retrieval.

Fig. 5 :

 5 Fig. 5: Graph partioning with K = 16 results on several images from the Caltech-101 benchmark. GC: Graph Cuts, NC: Normalized Cuts.

TABLE I :

 I Influence of parameter λ (equation 1) on image retrieval. For each value, the mean average precision is given. For this experiment, the patch size is set to n = 5.

		Mean Average Precision for the IMPK approach with 3 levels
		λ = 0	5	20	100
	NC	0.386079	0.389737	0.386185	0.380554
	GC	0.373073	0.382115	0.377872	0.375660

TABLE II :

 II Influence of parameter n (section III-B) on image retrieval. For each value, the mean average precision is given. For this experiment, λ = 5.

		Mean Average Precision for the IMPK approach with 3 levels
		n = 3	5	7	9
	NC	0.380546	0.388987	0.389456	0.388229
	GC	0.379723	0.382115	0.377215	0.375741

1 http://ai.stanford.edu/ bangpeng/ppmi.html

TABLE III :

 III Retrieval performance (mAP and standard deviation) on Caltech-101 dataset.

		Level 0 (BoW)	Level 1	Level 2	Pyramid
	NC	0.117 ±0.18	0.107 ±0.16	0.125 ±0.19	0.129 ±0.20
	GC	0.117 ±0.18	0.099 ±0.15	0.114 ±0.18	0.122 ±0.18
	SPM	0.117 ±0.18	0.144 ±0.21	0.162 ±0.23	0.157 ±0.22

TABLE IV :

 IV Retrieval performance (mAP and standard deviation) on PPMI dataset.

		Level 0 (BoW)	Level 1	Level 2	Pyramid
	NC	0.109 ±0.031	0.104 ±0.022	0.106 ±0.025	0.108 ±0.027
	GC	0.109 ±0.031	0.104 ±0.022	0.106 ±0.024	0.108 ±0.027
	SPM	0.109 ±0.031	0.115 ±0.036	0.121 ±0.043	0.118 ±0.041

http://trec.nist.gov/

ACKNOWLEDGMENT

This work was partially supported by CNRS (Centre national de la recherche scientifique) & Region of Aquitaine Grant.


A Bilbao

I Azpeitia

S Luyckx

N Loix

J Muñoa

Low Frequency Chatter Suppression using an Inertial Actuator

Keywords: Active damping, Chatter, Inertial actuator

It is well known that inertial actuators can remove regenerative chatter in many milling processes. They are useful when active control techniques are needed and the use of passive dampers is not enough. Moreover, complex control laws can easily be implemented leading to large improvements. However, the main limitation of inertial actuators is that they can only be used in a particular range of frequencies. This limitation is due to a low frequency suspension mode, which only allows the actuator to work at larger frequencies than this suspension mode. This is an important drawback in many big milling machines where chatter appears close to this low frequency. This paper deals with this problem and shows a possible solution to avoid this limitation from a control point of view. Some real milling tests will be presented showing the usefulness of the technique.

INTRODUCTION

Nowadays, self-excited vibrations or regenerative chatter is one of the main limitations of milling processes. Machine tool chatter vibrations avoid obtaining the required accuracy in workpiece, reducing the lifetime of the tool and the mechanical components of the machine.

Probably one of the most used solutions to increase the stability margin is by adding a tuned passive vibration absorber to the structure [START_REF] Sims | Particle dampers for workpiece chatter mitigation[END_REF][START_REF] Muñoa | Optimization of multiple tuned mass dampers to suppress machine tool chatter[END_REF]. This technique is not new and it has been widely used since it was proposed in [START_REF] Hartog | Mechanical vibrations[END_REF]. However, a passive absorber is not feasible in many machining processes where the dynamics of the system change during the process and then an active damper is needed [START_REF] Chung | Active damping of structural modes in high speed machine tools[END_REF][START_REF] Preumont | Vibration Control of active structures. An Introduction[END_REF][START_REF] Brecher | Electrohydraulic active damping system[END_REF][START_REF] Ganguli | Active damping of chatter in machine tools -demonstration with a hardware-in-the-loop siulator[END_REF]. Moreover, it is common to observe a non linear behaviour in passive absorbers making more difficult an accurate tuning of it.

Another possibility is to include an inertial actuator attached to structure to be controlled. An inertial actuator can easily be seen as a device producing a reaction force on the host structure according to a designed control law [START_REF] Huyanan | Vibration control strategies for proof-mass actuators[END_REF][START_REF] Sabalza | Incremento de la capacidad productiva de las fresadoras mediante la reducción del efecto del cahtter utilizando actuadores inerciales[END_REF]. For this reason, inertial actuators have been used to reduce vibrations in machining operations. The active damping is achieved by implementing an external control law. The vibration is measured on the structure and then the designed control law generates a control signal to the actuator. The main drawback of inertial actuators is that they usually possess a suspension mode at low frequencies. This means that actuators only behave as a linear force generator beyond its natural frequency. This fact could be a problem in many big machines where main modes appear at the same low frequency range. Depending on the materials used in their construction, a non-linear behaviour could be seen in many actuators as well.

This paper deals with this problem and shows a possible solution to this limitation. We are proposing an easy solution from a control point of view. Since the control law is implemented by software there is no need to change to a more appropriate actuator with a larger linear range. Many other solutions could be found proposed by other authors. In [START_REF] Paulisch | Active Vibration Control Using an Inertial Actuator with Internal Damping[END_REF] an internal damping is proposed to increase the stability of the controller. Similar results can be see if force feedback is implemented [START_REF] Bernassi | Active vibration isolation using an inertial actuator with local force feedback control[END_REF]. The main idea of this type of controllers consists of implementing a second control law dealing with the dynamics of the actuator. However, this is not always possible since some of the signals needed for the internal feedback could not be available for measuring. The use of a Kalman filter or a similar state-space observer could make an estimation of this signal available, but this is not always an easy solution. Another possibility is to implement more complex control laws. A Virtual Passive Absorber (VPA) controller is a good option since it does not excite the dynamics of the actuator [START_REF] Preumont | Vibration Control of active structures. An Introduction[END_REF][START_REF] Sims | Vibration absorbers for chatter suppression: A new Analytical tuning methodology[END_REF][START_REF] Wu | Adaptive control using a vibration absorber controller[END_REF]], but it is not very useful if the vibration chatter is close to the low frequency dynamics of the actuator. Another possible controller consists of implementing a displacement feedback in parallel with the main control law. A displacement feedback is able to affect the stiffness of the mode to be controlled and change its location into the linear range of the actuator. Then main control law can be applied.

In this paper we are including a compensation filter in the control-loop. This kind of filters can increase the working range of the actuator and avoid the suspension mode limitation. The objective of this kind of filters is to modify the dynamics of the actuator seen by the controller such that it does not interact with the structure dynamics anymore. This approach leads to good results. However, it relies on a good understanding of the actuator dynamics and, obviously, there are always non-modelled dynamics in the model of the actuator. The more accurate the model is, a better closed-loop performance is obtained.

Then a second drawback should be taken into account. Some inertial actuators possess a non-linear behaviour. This means that it is not possible to find a unique linear model describing the behaviour of the actuator for all the possible operation points. A set of different operation points is considered. Then a different transfer function is calculated for each operation point. Each of this transfer functions is calculated via linearisation of the behaviour of the actuator around the operation point. Then a compensation filter is associated to each transfer function. Finally, once the dynamics of the actuator are moved to lower frequencies, then a controller can be implemented to remove the undesired vibration. In this paper a direct velocity feedback has been used. This kind of control law is widely used on inertial actuators due to its simplicity and leads to good results. There is plenty of information about this control explaining its advantages and disadvantages [START_REF] Preumont | Vibration Control of active structures. An Introduction[END_REF][START_REF] Ehmann | Comparison of control strategies for active vibration control of flexible structures[END_REF]. This paper is organized as follows. First, the concept and a general definition of an inertial actuator is explained. Then the actuator used for the work is presented. The next section describes the control strategy, including the choice of the compensation filter and the direct velocity feedback control law. Some of the performed experiments will be presented showing the usefulness of the solution proposed in the paper. These tests include some real milling tests where the chatter vibration is removed by the actuator. Finally, conclusions and future work will end the paper.

INERTIAL ACTUATOR

An inertial actuator can be described as a reaction mass ma supported on a spring ka and a damper ca attached to a base (Figure 1). According to an applied voltage Vin the reaction mass is excited and it induces a force Fa on the supporting base. The dynamics of an inertial actuator can be modelled as the following transfer function

F a  s  V in  s  = G •m a s 2 m a s 2 c a s k a =g a s 2 s 2  2 a  a s a 2 (1)
where ωa is the natural frequency, ζa the damping ratio and ga is the gain of the actuator. The actuator behaves as an ideal linear force generator beyond a determined frequency ωp (Figure 2). Another upper limitation in the bandwidth is included by the electromagnetic circuit. This means that the actuator can only be used as a force generator in this linear range. In this paper the actuator used to perform the practical applications is a model of Micromega Dynamics (http://www.micromega-dynamics.com) (Figure 3). This actuator is capable of applying an up to 1000 N force on the host structure. Moreover, it is designed to work in two perpendicular axes. This is an important characteristic about this actuator, since most of them are designed to work in only one direction. In order to identify the behaviour of the actuator, it is mounted on a Kistler dynamometer plate and excited by a chirp voltage. Different amplitudes of the chirp signal are considered and the obtained frequency response functions (FRF) show a non-linear behaviour of the actuator (Figure 4). The natural frequency of the actuator depends on the amplitude of the chirp signal. Larger values of this amplitude lead to lower frequencies of the suspension mode. Most of the times this is not important if the vibration to be damped is in the linear range of the actuator. However, this non-linearity could affect the performance of the controller if the vibration to be remove is close to this suspension mode.

Direct Velocity Feedback

In this paper a direct velocity feedback is implemented to introduce some damping in the process. The transfer function of the controller can be written as

V c  s=-g v Ẋ  s=-g v • s • X  s (2)
where Vc(s) is the applied control input to the inertial actuator, s•X(s) is the velocity on the structure and

v g R + ∈
is the control gain. Considering that the actuator is working in the linear range, then the transfer function between the applied force on the structure and the control voltage input can be written as

F a  s =g a • V c  s =-g a • g v • s • X  s (3)
If the structure is modelled as a linear time invariant SISO system, the characteristic equation can be calculated as follows

1 g v • g a • s •G s=0 (4)
where G(s) is the transfer function describing the structure to be controlled. Now considering a structure with a unique mode, equation ( 4) can be written as

1 g v • g a • s • g s 2 2 n  n s n 2 =0  s 2  s  g v • g a • g2  n  n  n 2 =0 (5)
From ( 5) the damping of the closed loop poles can be augmented increasing the value of the feedback gain gv. However, the effect of the inertial actuator should taken into account. The dynamics of the inertial actuator can be studied including its transfer function in previous characteristic equation

1 g v •s • g a s 2 s 2 2 a  a s a 2 G  s =0 (6)
In this case the low frequency dynamics of the actuator should be considered to guarantee closed-loop stability.

In Figure 5 a root locus analysis is performed for a unique mode structure including the actuator dynamics. It shows that the poles of the structure are damped increasing the feedback gain. However, the poles of the actuator become unstable for large values of gv. This means that for sufficiently large values of gv the close loop is unstable. The largest value of this gain that keeps the system stable could be calculated if a linear model of the structure was available. However, this is not always possible and the optimal value of the feedback gain gv is usually calculated experimentally. Initially it is chosen to be small and it is increased until the closed-loop system becomes unstable.

Figure 5: Root locus analysis including the dynamics of the actuator

Compensation Filter

When the vibration to be removed is close to the suspension mode of the actuator, then a controller based on just a direct velocity feedback could not be enough to achieve a good result. The optimum control is achieved when the applied force by the actuator possesses the opposite phase than the velocity of the structure. The phase of the voltage applied to the actuator is always close to the correct one since it is generated by the controller. If the vibration to be damped is in the linear range of the actuator, then the applied force has similar phase to the voltage input and the desired effect is achieved. However, this does not happen if the vibration is close to the suspension mode. Then there is not certainty that the applied force has the correct phase. Different solutions can be applied to avoid this limitation. One could deal with the structure and try to move the main mode to higher frequencies, where the actuator behaves as a linear force generator. A very simple strategy is to include a position feedback controller. The effect of this controller is to move the frequency of the mode by changing its stiffness. Then, once both modes are sufficiently separated, one can apply the velocity feedback and expect some damping. The disadvantage of this strategy is that some amount of the force that the actuator can apply is wasted moving the mode to another frequency and full capacity can not be used in the damping process. Another solution is to increment the linear range of the actuator moving its suspension frequency to lower frequencies. This can easily be done including a compensation filter in the control-loop. A compensation filter modifies the dynamics of the actuator seen by the controller such that it does not interact with the structure dynamics anymore. Of course it relies on a good understanding of the actuator dynamics and in practice there is always a difference between the real behaviour and the linear model used to define the filter. This means that a perfect actuator dynamics is never achieved. However, most of the times this technique leads to good results and there is always some improvement.

One could think that it is not always easy to find a linear model describing the behaviour of the actuator. In previous section, it was shown that the actuator presents a non-linear behaviour. Then it is not possible to find a linear model describing its behaviour for any condition. In this paper, we are proposing to linearise the behaviour around the operation point. This technique requires an a priori knowledge of the process. However, this is a feasible condition for many milling processes, where it is common to know the value of the different signals of the process. Then a set of compensation filters are considered. Each one calculated from a different operation point in such a way that the dynamics of the actuator seen by the controller is good enough for the control objectives.

EXPERIMENTAL RESULTS

Impact Tests

In this section the experimental results are presented. The different tests are performed using a Soraluce milling machine The actuator should be mounted as much as closer to the tool, where the cutting process takes place. The actuator can not be implemented in the spindle head due to its size, so it is attached under the ram of the machine (Figure 6). Moreover, it was designed to compensate the effect of gravity on the moving mass when it is mounted in that position. Then the part of the ram closer to the spindle head, where the displacement is larger, is the best option.

Figure 6: The actuator mounted under the ram First a dynamical analysis is performed hitting the structure with an impact hammer. Then the response on the structure is measured with an accelerometer. The direct velocity control law is used. Different outputs for the ram are considered and both horizontal and vertical axes are tested. Figure 7 shows the frequency response functions (FRF) for different values of the feedback gain. The output of the ram is 1100 mm and the horizontal axis is first excited. The main mode of the structure can be damped by increasing the feedback gain. On the other hand, both actuator mode (22Hz) and machine suspension (15Hz) are easily excited. This effect limits the effect of the controller and prevents achieving a good performance. Note that gv is a constant that multiplies the velocity measured on the structure and the result is a voltage signal directly applied to the actuator, so its units are V•m -1 •s.

Figure 7: FRF for different values of the feedback gain In order to improve this result a compensation filter (C1) is included in the control loop. The filter is designed to locate the suspension mode of the actuator at 10 Hz. In this case, since an impact is used, the linearisation of the behaviour of the actuator for low levels of the excitation has been used. Then different values of the feedback gain are considered and the FRFs are displayed in Figure 8. The mode of the actuator is not excited, but it appears an oscillation at 15 Hz. This is due to the fact, while the new actuator dynamics is now able to damp the main machine mode, it is not yet good enough to tackle the other mode of the machine, which is close to the new dynamics of the actuator.

Figure 8: FRFs for different feedback gains using the compensator filter C1 It is of course possible to improve the compensation filter in order to be able to address also the machine mode at 15 Hz. The compensation filter is chosen to locate now the suspension mode of the actuator at 5 Hz. Figure 9 shows the FRF for the new compensation filter (C2). Now the oscillation appears at lower frequencies, but with this compensator larger values of the feedback gain can be used before the instability appears. Moreover, the mode at 15 Hz is not excited. It is also possible to consider that there is no need to damp the machine mode at 15 Hz since it has little influence on the chatter vibration. Then the original compensation filter can be used if a "Notch filter" is added at 15Hz to prevent the control to address this machine mode.

Milling tests

Finally, some milling test are performed in the same milling machine (Figure 11). The cutting conditions are displayed in Table 1, while the cutting tool is described in Table 2. Note that the experiments will be performed for two different values of Fz. The compensation filter is chosen considering that the actuator needs to apply a large force to cancel the chatter vibrations on the ram. The results without the compensation filter are not shown since the actuator was easily excited and large values of the feedback gain could not be used. First Fz = 0.2 mm/z is chosen. Figure 12 shows the stability conditions for different values of the ram output and the depth of cut. The real values of the depth of cut are not given due to confidentiality concerns and all the values are compared to the minimum value of the depth of cut leading to an unstable process (100%). For lower values of ap, the process is always stable and the actuator is not needed to avoid the appearance of chatter. However, larger values than 100-125% lead to unstable conditions. Then the controller is connected and the stability margin is doubled for almost all the ram outputs. This means that the productivity of the machine is increased as well. The MRR (material removal rate) has been doubled. For very large values of ap the actuator is not able to remove completely the chatter vibration. However, there is always some improvement and the vibration is always reduced when the controller is connected. The cutting conditions are changed reducing the feed per tooth to 0.1 mm/z. In this case the controller actuator is able to remove the chatter vibrations for the maximum depth of cut and all the ram outputs (Figure 13). This means that vibrations do not limit the productivity and any condition can be used. This is always the best result when active control is applied in milling processes. One could ask why the stability margin could not be increased until the maximum value when Fz = 0.2 mm/z. The most common limitations are the dynamics of the actuator or a saturation of the force applied by the actuator. In this case, the limitation comes from the fixturing table used in the process (Figure 11). A new test is performed with ap larger than 200%. According to Figure 12 this condition leads to an unstable process. In Figure 14 the acceleration measured during the test is displayed. First the actuator is switched off resulting to an unstable process. After a few seconds the actuator is connected with gv = 2200 V•m -1 •s and the vibration on the ram is reduced. A frequency analysis is performed in order to obtain more information about the process (Figure 16). When no controller is applied then a chatter vibration appears at 34.5 Hz. Clearly this is due to the structure of the milling machine since its main mode is close to this frequency. Then when the controller is connected, this vibration is removed, but it appears another at 46 Hz, which results to be close to the main mode of the fixturing table (48Hz). This effect is more evident if the acceleration measured on the fixture table is analysed. Figure 15 shows the acceleration of it during the same test. One could see how the signal is increased when the controller is applied. Moreover, a frequency analysis of the velocity (Figure 17) leads to the same conclusion. The main conclusion from this result is that the machine is no longer the most important limitation; the fixturing table is now the main limitation. Different solutions can be implemented to deal with this drawback. Changing the fixturing table is the most evident one, but similar control techniques can be applied as well. A passive absorber or a second actuator could be attached to the fixturing table.

CONCLUSIONS AND FUTURE WORK

In this paper an active control using an inertial actuator for low frequency chatter suppression in milling processes has been presented. The appearance of low frequency chatter limits the performance of most inertial actuator due to the suspension mode. An easy solution from a control point of view has been used leading to a lower sensitivity to chatter instabilities. The inclusion of a compensation filter in the control loop increases the linear range of the actuator moving the suspension mode to lower frequencies. Moreover, the non-linear behaviour of the actuator has been overcome by considering a set of possible linearisation of the actuator around different operating points.

Results show the necessity of an a priori knowledge of the conditions for an appropriate choice of the compensation filter. The performed milling tests show that the productivity and the stability margin can be increased by the actuator even if the chatter frequency is close to the suspension mode of the actuator. The machine was found to be no more the main limitation for the appearance of chatter.

On the other hand, the fixturing table used for the tests introduces an important limitation. The future work will be focused on the solution to this drawback. The main idea consists of adding an additional absorber to the fixturing table. This absorber could be a passive absorber or a second inertial actuator attached to the fixturing table.

Figure 1 :

 1 Figure 1: Inertial actuator

Figure 2 :

 2 Figure 2: Bode plot of an ideal inertial actuator

Figure 3 :Figure 4 :

 34 Figure 3: Inertial actuator from Micromega Dynamics

Figure 9 :Figure 10 :

 910 Figure 9: FRFs for different feedback gains using the compensator filter C2

Figure 11 :

 11 Figure 11: Set up for milling tests

Figure 12 :Figure 13 :

 1213 Figure 12: stability conditions for different ram outputs and depths of cut when Fz = 0.2 mm/z

Figure 14 :Figure 16 :Figure 17 :

 141617 Figure 14: Acceleration measured on the ram with ap larger than 200%

Table 1 :

 1 Cutting Conditions

	Cutting tool
	Reference	SANDVIK R245-125Q40-12M
	D (diameter)	125 mm
	Z (inserts)	8
	k (lead angle)	45 degrees

Table 2 :

 2 Cutting Tool

ACKNOWLEDGEMENTS

The authors are very grateful to the European Commission by its support through Chameleon project CP-FP7-NMP2-SL-2008-213319-2 and to the Department of Industry, Innovation, Commerce and Tourism of the Basque Government by its support through contract IEI09-035 'Hiperion-class heavy-duty machines' of ETORGAI program.