Sughosh Deshpande

Maria Clara

Coimbra Gonçalves

Anna Carla Araujo
email: araujo@insa-toulouse.fr

Pierre Lagarrigue

Yann Landon

Specific force map for smart machining applications with rotating tools

Keywords: Material identification, monitoring, helical milling, axial drilling, cutting forces, specific coefficients 1

Industry 4.0 is the need of the hour in current global market scenario and all the processes are moving towards automation and smart manufacturing. In machining, smart techniques implementation depends on developing a database for decision-making, which is the case for stack drilling in aerospace industry. In this application, choosing one optimal condition for several materials is a challenge due to their different machinability. Hence, material identification techniques are suitable approaches for adapting the cutting parameters in real time, which improves tool life, hole quality and productivity. In that regard, the goal of the present paper is to create a specific force data map for axial drilling and circular milling processes based on its experimental force and power measurements. To do that, experiments were separately carried out on Titanium and Aluminium workpieces in a range of cutting speed and feed conditions. The results show that specific cutting and feed forces for each material can be identified on distinct regions of the map, without thresholds overlapping. Given that, these maps can be used as a signature to distinguish two metallic materials in real time machining. In this case, the specific data points at the interface layers may offer advantage to accurately identify tool position unlike monitoring gradient of feed forces while drilling stacked materials. Therefore, smart machining techniques seeking cutting parameters optimization can be implemented for a particular material.

Introduction

Multi stack drilling is one of the major operations in aircraft assemblies and it demands lots of time and resources to achieve high-quality holes.Usually, Titanium and Aluminium stacked materials are drilled before aircraft assembly and this poses a major challenge in machining due to the different materials' machinability. Some of the drilling challenges include poor finishing, burr formation and rapid tool wear [START_REF] Zhu | Evaluation of novel tool geometries in dry drilling aluminium 2024-T351/titanium Ti6Al4V stack[END_REF]. These hole quality aspects, as exit burrs, can have a significant influence during assembly and affect productivity [START_REF] Gillespie | Burrs produced by drilling[END_REF]. Hence, as thousands of holes have to be drilled for assembling an airplane, techniques considering different materials' machinability in a stack are crucial for performance maximization.

In this aspect, smart machining is a solution that refers to real time adaptation of cutting parameters for process optimization based in massive data from smart factories, specially using machine tool acquisition [START_REF] Kong | A study on a general cyber machine tools monitoring system in smart factories[END_REF]. Previous research have presented prediction models and data monitoring techniques which can be helpful in smart machining applications. Wenkler et al. [START_REF] Wenkler | Intelligent characteristic value determination for cutting processes based on machine learning[END_REF] developed a way of predicting specific cutting force by Artificial Neutral Network which can be implemented in process planning or smart manufacturing. Hegab et al. [START_REF] Hegab | A smart tool wear prediction model in drilling of woven composites[END_REF] predicts tool wear by using Machine Learning models for a wide range of cutting conditions in drilling. Furthermore, a smart machining system is developed by monitoring spindle power and torque in order to optimize feed rates in milling by Park et al. [START_REF] Park | Development of smart machining system for optimizing feedrates to minimize machining time[END_REF].

One interest of applying smart machining during stack drilling is material identification. With this purpose, tool position detection inside the stacked layers is a key parameter so it is possible to adapt cutting parameters and machining strategies. Regarding this aspect, Pardo et al. [START_REF] Pardo | Assessment of decision-making algorithms for adaptive drilling of aerospace stacks[END_REF] used cutting forces signals as input for three different decision-making algorithms to identify the tool position in stacked materials. Their objective was to compare the effectiveness of different approaches, which could be used for material identification. Jallageas et al. [START_REF] Jallageas | Self-adjusting cutting parameter technique for drilling multi-stacked material[END_REF] developed a methodology to material identification in drilling multi-stacked material using thrust force signals. Their strategy used both: thresholds and gradient of forces to differentiate metallic materials and spectral analysis in order to distinguish aluminum and carbon fiber. The approach does not use any input process parameter, it uses pattern recognition for a specific drilling case. Material identification in stack drilling for a wide range of parameters is still a challenge that requires more research and the development. However, implementation of smart machining strategies in drilling alone may not completely optimize the process, so additionally other alternatives such as helical milling, also known as orbital drilling, can be applied.

Orbital drilling shows advantages regarding the process flexibility and ability to produce high-quality holes on different materials [START_REF] Bruno | A review of helical milling process[END_REF][START_REF] Morales | Study of the Influence of Helical Milling Parameters on the Quality of Holes in the UNS R56400 Alloy[END_REF]. Sun et al. [START_REF] Sun | Hole-making processes and their impacts on the microstructure and fatigue response of aircraft alloys[END_REF] data shows that the fatigue life of holes produced by helical milling showed better life compared to conventional drilling. Although, cutting mechanistic is more complex in this process and so is the cutting models. Wu et al. [START_REF] Wu | Cutting force prediction for circular end milling process[END_REF] discussed force models for circular milling taking into account the uncut chip thickness and feed rates variation at the corners. The cutting coefficients calculation in circular model can be further extended for helical milling adding vertical feed.

In both hole making processes of aerospace materials: axial drilling and helical milling, the specific cutting coefficients data can be used to identify the material being machined. They can act as material signatures to adapt proper cutting parameters in real time. These data could directly improve the process productivity and hole quality if appropriate cutting parameters are applied for the specific material.

Hence, the goal of this paper is to develop a data map of specific force coefficients related to each material for both process: axial drilling and circular milling processes, not developed in the previous papers, for further application in helical milling. The materials used for the experiments were aluminium 2017A and titanium Ti6Al4V, being subjected to a common set of cutting conditions. In the following sections, the force model and cutting coefficients are presented, as well the experimental set-up. The cutting forces and power results are presented and discussed, and lastly, the identification map with cutting coefficients for both processes is built.

Force model and specific force coefficients

Mechanistic force models can provide quantitative cutting force predictions based on the uncut chip thickness principle. Once experimental tests are made, the cutting force models enable estimating the specific cutting coefficients for a different set of variables, as describe by numerous papers [START_REF] Armarego | Computerized Predictive Cutting Models for Forces in End-Milling Including Eccentricity Effects[END_REF][START_REF] Altintas | Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design[END_REF].

The forces for a general orthogonal cutting can be described by a mechanistic force model considering only cutting action of the cutting edge neglecting ploughing and chisel edge effects. The machining force ⃗ F m is calculated considering the small finite elements of the cutting flute in the cutting edge referential frame decomposed: tangential (dF c ), radial (dF r ) and axial (dF z ) components:

⃗ dF m =   dF c dF r dF z   . ( 1
)
The local force in point P of the cutting edge is a function of the uncut chip thickness h and the specific force coefficients: K c , K r and K z . dF c = K c (P ).h(P )db dF r = K r (P ).h(P )db dF z = K z (P ).h(P )db where h(P )db is the elementary uncut chip load calculated on the reference point P in db elemental length of the cutting edge.

For each rotating tool machining process there are some assumptions in order to identify a global value for specific forces coefficient: in drilling, milling and circular milling, as follows.

Drilling

In drilling, local specific force coefficients can change along the cutting edge, as cutting speed varies from zero to the maximum cutting speed V c on the tool nominal radius D/2. Also, cutting angles can change along the edge. Although, for a fixed tool geometry, it is possible to calculate an average specific force coefficient along the edge. It can be calculated using the total chip load per cutting edge

A c = ∫ hdb = f z (D/2) = f 2 (D/2
) (where f is the feed per revolution for a 2 flutes tool and f z is the feed per tooth).

In drilling, it is not possible to experimentally measure the cutting force for a single cutting edge F c (N), because both edges are in contact with the workpiece and they have opposite directions. The average specific cutting force K c (N/mm 2 ) for a tool-workpiece pair in drilling is calculated using the experimental cutting power P c , contribution of both cutting edges [START_REF] Wenkler | Intelligent characteristic value determination for cutting processes based on machine learning[END_REF]:

K c = F c A c = 120 P c V c f z D (2)
On other hand, it is possible to measure the feed force F z and to directly calculate the feed force F f for each cutting edge (F z = 2F f ). The specific feed force K f (MPa) is calculated as:

K f = F f A c = F f f z (D/2) = F z f z D (3)

Milling

In milling, in order to determine the cutting coefficients of the tool-workpiece pair, the local cutting force dF c (t) and the radial force dF r (t) are analyzed considering the variation of the uncut chip thickness along the tool rotation.

Martellotti equation [START_REF] Araujo | A model for thread milling cutting forces[END_REF] described h(θ 2 ) as a function of the tool rotation angle θ 2 . The forces in axial direction (dF z ) can be neglected, if the helix angle is small. When summing the contributions of each cutting part db, dF c (t) and dF r (t), the directions are not the same and a change to a fixed referential is necessary. From there, dF x (t) and dF y (t) can be calculated and the machining tool for each cutting edge in contact to the workpiece at the time t can be calculated.

As in drilling, some assumptions should be made for the identification of specific cutting force and specific radial force. It is considered that the maximum force is achieved in the position of θ 2 when the average uncut chip thickness along the cutting edge is higher. The reference frame is now fixed in the rotating tool θ 2 and not in each element. Then, the approximation done lead to:

K c = F c (θ 2 ) A c (θ 2 ) -→ K c = F c (max) A c (max) ≈ F c (max) h ′ (max).a p (4) K r = F r (θ 2 ) A c (θ 2 ) -→ K r = F r (max) A c (max) ≈ F r (max) h ′ (max).a p (5)
Considering that the measurement is taken using a fixed dynamometer, the acquisition is done in a XY fixed referential frame. In linear milling, when feed direction is aligned with X direction :

[

F X F Y ] = [ Sin(-θ 2 ) -Cos(-θ 2 ) Cos(-θ 2 ) Sin(-θ 2 )
] .

[

F c (θ 2 ) F r (θ 2 ) ] ( 6
) X Y F t F r F t F a F r F tan F rad

Circular Milling

Unlike in linear milling, circular milling involves movement of tool around a centered fixed point in clockwise or anticlockwise direction to enlarge a existing hole diameter. Figure 1(b) shows the circular milling tool trajectory (at θ 1 = 0) and the associated geometric parameters. A tool of diameter D t with center C is moving around the hole center O. D 1 refers to pre hole diameter and D is the final diameter to be achieved after circular milling. R tt is the tool trajectory radius with reference to hole center O [START_REF] Araujo | Analytical and experimental investigations on thread milling forces in titanium alloy[END_REF].

In circular milling, cutting forces at the cutting edge have different referential frames that have to be analyzed in order to determine the cutting coefficients of the tool-workpiece pair. The cutting forces exerted by the tool on the workpiece can be seen in Figure 1(a). Feed motion produces F tan and F rad forces and the cutting action gives rise to tangential (F t ), radial (F r ) and axial forces (F z ) at the cutting edge. It can be noted that the direction of feed forces in circular milling changes continuously unlike in linear milling where feed direction and feed force direction remain unchanged. In circular milling context, forces in axial direction (F z ) can be neglected.

The forces measured during experiments by the dynamometer are in XY referential frame: F X and F Y as shown in Figure 1(a). The tool rotates around the hole center O in clockwise direction given by angular position θ 1 and around its own axis given by θ 2 , also in clockwise direction. The forces F X and F Y can be transformed using a rotation matrix (Eq. 7) to deduce F tan and F rad at the tool center. F tan and F rad are significant to understand tool deflection effects on hole quality. To calculate F t and F r at the cutting edge, the rotation matrix as a function of θ 2 and the dynamometer referential frame components (F X and F Y ) (Eq. 8) can be applied. The specific cutting coefficients K c and K r (N/mm 2 ) are identified by analyzing F t and F r maximum values for a particular set of tool revolutions using Eq.4 and Eq.5.

[

F rad F tan ] = [ Cos(-θ 1 ) -Sin(-θ 1 ) Sin(-θ 1 ) -Cos(-θ 1 )
] .

[

F X F Y ] (7)
[ F t (θ 2 ) F r (θ 2 ) ] = [ -Sin(-θ 2 ) -Cos(-θ 2 ) -Cos(-θ 2 ) -Sin(-θ 2 )
] .

[

F X F Y ] (8)
3 Materials and methods

In this section, drilling and circular milling experiments are presented in Aluminium (2017A) alloy and Titanium (Ti6Al4V) alloy workpieces measuring cutting forces and cutting power. The experiments were carried out for a specific range of cutting conditions, following the experimental setup and the design of experiments described below.

Experimental setup

All drilling and milling experiments were carried out on a CNC milling center DMU85-DMG mono block machine. Flood water-based through coolant was used for drilling and circular milling was carried out in dry condition. The tools used for the experiments were: carbide drills from Klenk having diameter 11 mm, one of the most commonly used in aerospace stack drilling applications, and carbide end mills of 8 mm diameter from Fraisa.

The workpieces were fixed on a 9257B Kistler dynamometer, as shown in Analogical data of force and power were converted to digital using a 9201 National Instruments acquisition module with 10 kHz acquisition rate. The measured forces were filtered using a low pass Butter-worth filter with cutoff frequency of 700 Hz (appropriated to the maximum spindle speed).

Design of Experiments

Table 1 summarizes the cutting conditions on Aluminium and Titanium along with tool details for drilling and circular milling operations. The range of cutting parameters is wider than usual because the goal was to compare the responses in both materials and processes for a wide range of cutting conditions. The lowest feed per tooth allowed for a monoblock carbide tool is included (0.02mm/tooth) in the parameters range.Drilling depth was set as 15mm in case of axial drilling for both materials. Circular milling was carried out on the drilled hole with an axial depth of 1mm and a radial width of cut of 1.9mm. In order to have a cutting continuity coefficient c less than 1, in other words to have only one cutting flute in contact with the workpiece at any point of time [START_REF] Araujo | Modeling thread milling forces in mini-hole in dental metallic materials[END_REF], the chosen final milled hole diameter is 14.8mm.

Experimental results

The following subsections present the results of feed force and power for drilling along with cutting forces for circular milling. The graphs are presented in terms of length of the hole for drilling for various cutting speeds and feed rates. Circular milling forces are presented in terms of tool angular position θ 1 (around O) which is taken as a reference for all the feed rates and cutting speed values tested. It should be noted that the trajectory rate (θ 1 /time) changes based on variation in feed rate values.

Feed force and power in drilling

Regarding the cutting force, Figure 3 shows the results obtained per depth of the hole for Aluminium and Titanium respectively. Each figure contains results of different cutting speeds obtained at f z =0.1 mm/th. Varying cutting speeds do not greatly differ the feed force values as seen in the Figure 3. Also, the region where the tip of the tool starts entering the workpiece is clearly observed from the increasing force slope, then the cut is stabilized. The complete graphical results from the experiment for other feed values are referred in the earlier article from Gonçalves et al. [START_REF] Clara | Smart drilling: material identification using specific force map[END_REF].

Cutting power results are shown in Figure 4 for Aluminium and Titanium at f z =0.1 mm/th. For Titanium, there is an increase in power when the tool penetrates the workpiece unlike Aluminium which remains almost stable during drilling. This can be attributed to the elastic recovery of material, since no tool conicity was observed. An important fact that can be observed comparing both graphs is: the results are the same during the first 2 mm penetration for both materials, when only the chisel point is in contact with the workpiece. In consequence, identifying the materials during this phase of drilling, and also in the transition of materials, is challenging in stack drilling. This was also observed previously by Pardo et al. [START_REF] Pardo | Assessment of decision-making algorithms for adaptive drilling of aerospace stacks[END_REF].

Forces in Circular Milling

The cutting forces F X and F Y obtained with the dynamometer during circular milling can be seen in the Figure 5. The graph shows the variation of forces Figure 6 shows the variation of transformed forces at the cutting edge (F t and F r ) in a time window of two tool revolutions θ 2 = [0 -720 o ] in Aluminium and Titanium at f z = 0.1 mm/th and V c = 40 m/min. Eight peaks can be observed in the transformed forces for a 4 flutes end mill. The cutting forces were transformed as previously indicated in Eq. 8.

Experimental specific force coefficients

After acquiring forces and power data for the analyzed processes, the specific cutting force coefficients were estimated. They are presented in the following subsections for the drilling and milling operations, followed by the targeted identification map. It is important to claim that the circular milling experiments were done after the analysis of drilling experiments. As the aim of this paper is the material identification based on the data, the range of feed per tooth was 10 Open Access Version of the articles Submitted to Part B: Journal of Eng. Manufacture in October-2021 narrowed to 0.10 and 0.18 mm/th, where the values of forces were closer. Also, low feed rate values as 0.02 mm/tooth are very rarely used in the industry for these materials.

Specific force coefficients in drilling

The K c and K f values for different feed rates are shown in Figure 7(a) and 7(b), where comparisons between both materials are made. The behavior of results are similar, considering that the values for Titanium are higher due to its poorest machinability. K c and K f values are calculated as indicated in the Eq.2 and 3 previously. The data includes values calculated at several tool angle positions θ 1 for each specific feed rate and cutting velocity with the objective of having a representative value with larger spectrum of experimental variability and noise.

It can be noted in Figure 8 that the specific coefficient values for Titanium are considerably higher compared to Aluminium, as observed for drilling.

Identification map

In this article we propose a clear distinction between both materials and force coefficient values identified in a single graphic, so called identification map. Even though the analysis of force coefficients per feed shows distinct results for different materials, identification can be mistaken through these results. This is because of the possibility of having similar coefficient values for both the materials at different feed rates.

Figures 9(a) and 9(b) presents the identification map using K c per K f for drilling and K c per K r for circular milling. It should be noted that in drilling K f is derived from the feed force and in circular milling K r is derived from the radial cutting force. The areas identified as Aluminium and Titanium are clusters with a simple method using minimum and maximum points on the map, without thresholds overlap. Although, the methodology could be improved using other sophisticated clustering methods. This clear distinction increases the material identification precision, which is the first step of self-adjustment techniques (changing spindle speed or feed per tooth).

In extreme situations, the acquired machining data can be incorporated to the commercial software in to alarm or stop the tool rotation in extreme situations, as it is already the case of power [START_REF] Köttner | Process monitoring using machine learning for semi-automatic drilling of rivet holes in the aerospace industry[END_REF].

Conclusions

The main focus of this article was to demonstrate the possibility of material identification during hole making operations of two different materials by analyzing specific force coefficients. In this regard, cutting force and power data were analyzed to identify specific force coefficients for Titanium and Aluminium. Our work proposes utilisation of data points of specific force coefficients for material identification while drilling stacked materials. This novel technique can help to adapt proper cutting parameters in real time. Such data maps with precise points of specific force coefficients can have advantages over monitoring of thrust forces for gradient change to identify materials or tool position in the The results can be summarized as follows:

• In both operations, specific force coefficients are identified in different regions of the map for Titanium and Aluminium due to the materials distinct machinability;

• During initial phase of axial drilling, when only chisel point is inside the material, the signals are similar comparing both materials and difficult to identify;

• Considering drilling results in Titanium, the elastic recovery of material leads to an increase in power when the tool penetrates the workpiece. This elastic recovery disturbs the identification, it could lead to wrong interpretation of increasing specific cutting force inside the hole, which is not the case.

• In circular milling, the transformed feed forces (F tan and F rad ) can be further analyzed to understand its relation with respect to hole quality and attributed to optimal cutting conditions for a tool-material pair.

• The variation of K r values in circular milling compared to K c values for a particular feed and velocity may be attributed to tool deflection/vibrations which is also reflected in radial forces. However, there is a distinct separation when the values are compared with both materials.

Hence, the gap between the cutting force coefficients for both materials makes it suitable for applying smart machining techniques using data monitoring in hole making processes of stacked materials made out of Aluminium and Titanium.

13

Open Access Version of the articles Submitted to Part B: Journal of Eng. Manufacture in October-2021

Figure 1 :

 1 Figure 1: Referential Frames in X and Y indicating cutting forces and tool trajectory in circular milling

Figure 2 :

 2 Figure 2: Experimental setup

Figure 3 :

 3 Figure 3: Experimental feed force for drilling Aluminium and Titanium

  Figure 2(a), using only its internal measuring region, connected to a 5070 Kistler m/min Vc = 40 m/min Vc = 20 m/min (b) fz = 0.10 (Titanium)

Figure 4 :

 4 Figure 4: Experimental cutting power for drilling Aluminium and Titanium

Figure 5 :

 5 Figure 5: Experimental forces F X and F Y during Circular Milling of Aluminium as an example (f z -0.1 mm/th and V c -40m/min)

  Fr for Titanium.

Figure 6 :

 6 Figure 6: Calculated F t and F r during 2 revolutions in circular milling of Aluminium and Titanium (window on θ 1 = [185 o -198 o ]) at f z = 0.1 mm/th and V c = 40 m/min

Figure 7 :

 7 Figure 7: Specific force coefficients (K c and K f ) for Aluminium and Titanium in drilling

Figure 8 (

 8 Figure8(a) and 8(b) shows K c and K r values against feed for Aluminium and Titanium. They were calculated using F t and F r values, respectively, for a specific set of tool revolutions.The data includes values calculated at several tool angle positions θ 1 for each specific feed rate and cutting velocity with the objective of having a representative value with larger spectrum of experimental variability and noise.It can be noted in Figure8that the specific coefficient values for Titanium are considerably higher compared to Aluminium, as observed for drilling.

Figure 8 :

 8 Figure 8: Specific force coefficients (K c and K r ) for Aluminium and Titanium in circular milling

  (a) Kc and K f (Drilling) (b) Kc and Kr (Milling)

Figure 9 :

 9 Figure 9: Identification Map for Drilling and Circular Milling in Aluminium and Titanium

Table 1

 1

					: Design of Experiments: Cutting conditions
	Operation						Workpiece	Cutting speed	Feed (mm/th)
										(m/min)
	Drilling							Aluminium 20 / 60 / 100	0.02 / 0.10 / 0.18
	( 11 mm, 2 flutes)				Titanium	20 / 40 / 60	0.02 / 0.06 / 0.10
	Milling							Aluminium 40	0.02 / 0.10 / 0.18
	( 8 mm, 4 flutes, 30°helix) Titanium	40	0.02 / 0.10 / 0.18
	0	2	4	6	8	10	12	14	16
				Length (mm)

Thus, as shown in this work, the specific force coefficients makes it possible to clearly identify different materials in real time, which is calculated by online monitoring of cutting forces and power. The calculation of instantaneous specific force coefficients for material identification and real time decision making to adapt cutting parameters will optimise the hole making process.


Qianyi Jiang
email: qianyi.jiang@eleves.ec-nantes.fr

Saïd Moussaoui
email: said.moussaoui@ec-nantes.fr

Jérôme Idier
email: jerome.idier@cnrs.fr

Guylaine Collewet
email: guylaine.collewet@irstea.fr

Mai Xu
email: maixu@buaa.edu.cn

Majorization-Minimization Algorithms for Maximum Likelihood Estimation of Magnetic Resonance Images

Keywords: Magnetic resonance imaging, Rician noise, maximum likelihood estimation, iterative optimization, majorizationminimization

This paper addresses maximum likelihood estimation of images corrupted by a Rician noise, with the aim to propose an efficient optimization method. The application example is the restoration of magnetic resonance images. Starting from the fact that the criterion to minimize is non-convex but unimodal, the main contribution of this work is to propose an optimization scheme based on the majorization-minimization framework after introducing a variable change allowing to get a strictly convex criterion. The resulting descent algorithm is compared to the classical MM descent algorithm and its performances are assessed using synthetic and real MR images. Finally, by combining these two MM algorithms, two optimization strategies are proposed to improve the numerical efficiency of the image restoration for any signal-to-noise ratio.

I. INTRODUCTION

Noise in Magnetic Resonance Imaging (MRI) is a major concern which may affect image interpretation and exploitation. In order to obtain images with high signal-to-noiseratio (SNR), long measurement duration is usually required. In practical use, it is necessary to reduce the acquisition time by associating some numerical noise reduction routines to the imaging system. Various kinds of denoising methods have been investigated (See for instance [START_REF] Mohan | A survey on the magnetic resonance image denoising methods[END_REF] which gives a thorough overview of existing methods) from which modelbased statistical estimation is becoming increasingly prevalent. In the case where the MR data are acquired in the complex domain (K-space), the noise in both real and imaginary parts is independent and identically distributed and it is usually supposed to be Gaussian [START_REF] Gudbjartsson | The Rician distribution of noisy MRI data[END_REF]. However, in most applications, the signal magnitude is processed since the phase is sensitive to many factors such as the magnetic field heterogeneity, the temperature of the sample or the motion of the analyzed object [START_REF] Peter B Roemer | The NMR phased array[END_REF]. Accordingly, the computation of magnitude raw images yields Rician distributed data which makes the denoising techniques based on the hypothesis of additive Gaussian noise suboptimal.

In the literature, several Rician denoising methods have been mentioned, including the so-called conventional approach, exploiting the properties of the second-order moment of the Rician distribution [START_REF] Mcgibney | An unbiased signal-to-noise ratio measure for magnetic resonance images[END_REF], the linear minimum mean square error (LMMSE) estimator [START_REF] Aja-Fernández | Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach[END_REF] and the signal reconstruction method based on maximum likelihood (ML) estimation [START_REF] Sijbers | Maximum-likelihood estimation of Rician distribution parameters[END_REF]. When the noisy data can be modeled accurately, the ML estimator is commonly applied since it tends to be consistent and asymptotically efficient under general conditions [START_REF] Ole | Parameter estimation from Rician-distributed data sets using a maximum likelihood estimator: Application to T1 and perfusion measurements[END_REF].

The negative log-likelihood (NLL) function for statistical estimation from Rician distribution is non-convex but unimodal. Although the conventional optimization methods can converge to the global minimum, their convergence efficiency depends severely on the initialization and the value to be estimated. In [START_REF] Getreuer | A variational model for the restoration of MR images corrupted by blur and rician noise[END_REF], Getreuer et al. proposed a maximum a posteriori (MAP) estimation with total variation prior and considered a somewhat complex convex approximation to overcome the non-convexity. Dong et al. [START_REF] Dong | A convex variational model for restoring blurred images with multiplicative noise[END_REF] established another convex Rician model by adding a quadratic penalty term into the previous model under some certain conditions to ensure the existence and the uniqueness of the solution. [START_REF] Solo | An EM algorithm for Rician fMRI activation detection[END_REF] and [START_REF] Zhu | Regression models for identifying noise sources in magnetic resonance images[END_REF] applied an Expectation-Maximization (EM) algorithm for Rician denoising estimation based on a convex tangent majorant function. Varadarajan [START_REF] Varadarajan | A majorize-minimize framework for Rician and non-central Chi MR images[END_REF] generalized this approach into the Majorization-Minimization (MM) framework for noncentral chi (NCC) family and demonstrated the computational efficiency compared with generic numerical optimization techniques. However, as presented in the following section, this algorithm results in a poor convergence speed for low values of signal-to-noise ratio (SNR) which is partly due to the nonconvexity of the NLL function. In this paper, we propose a novel MM optimization algorithm of the NLL function based on a transformed convex criterion that substantially accelerates the convergence speed at low SNR. By combining the conventional and our proposed MM approaches, we obtain two algorithms which improve the computational efficiency whatever the SNR level. The rest of the paper is organized as follows. Section II introduces the ML estimation in Rician noise. Section III overviews the principle of the existing non-convex model based MM algorithm. In Section IV, we introduce the novel MM algorithm and demonstrate the validation of our proposed majorant function. Section V compares the two algorithms and considers some combination strategies. Section VI illustrates the experimental results and Section VII presents some conclusions and perspectives of this work.

II. ML ESTIMATION FROM RICIAN DISTRIBUTION

Complex MRI signal is corrupted by an additive zeromean Gaussian noise. Thus, its magnitude follows a Rician distribution:

p (M i | A) = M i σ 2 e -M 2 i +A 2 2σ 2 I 0 AM i σ 2 u (M i ) (1)
where I 0 is the modified zeroth order Bessel function of the first kind, M i represents the ith measurement with the corresponding true signal amplitude A 1 , σ 2 is the variance of the Gaussian noise and u denotes the step function.

Under the assumption of independence of N measurements M i , the estimation of the pixel intensity A can be obtained by maximizing the likelihood function or equivalently minimizing the NLL function, given as:

L(A) = -log N i=1 p(M i |A) = N i=1 A 2 + M 2 i 2σ 2 - N i=1 log I 0 AM i σ 2 - N i=1 log M i σ 2 .
(

) 2
In the case where a single image is acquired (N = 1) hereafter for each amplitude A, there is only one measurement denoted M . In the sequel, the general case will be considered, and a discussion on the impact of the number of image accumulations will be given in section VI-A. As mentioned in [START_REF] Jan | Maximum likelihood estimation of signal amplitude and noise variance from MR data[END_REF], the objective function (2) has one or two stationary points where A = 0 is always a trivial solution.

III. MM DESCENT ON THE INITIAL CRITERION

The objective function cannot be minimized explicitly, therefore, an iterative optimization algorithm based on the MM framework can be applied [START_REF] David | A tutorial on MM algorithms[END_REF]. The principle of MM is to iteratively construct and minimize a surrogate cost function H L (•, A k ) known as tangent majorant instead of directly 1 The corresponding pixel index is omitted to simplify the notations. minimizing the initial objective function L(•). Each iteration k of the algorithm consists of:

Majorization: H L (A, A k ) L(A), ∀A 0 (3) Minimization: A k+1 = arg min A H L (A, A k ) (4)
The main feature of this descent method is that H L (•, A k ) admits a closed-form solution which is easier to compute than that of L(•).

The MM algorithm proposed in [START_REF] Varadarajan | A majorize-minimize framework for Rician and non-central Chi MR images[END_REF] takes into account that

f i (A) = -log I 0 c i A with c i = Mi σ 2 , is a concave function.
It is majorized at any point A k by an affine function:

f i (A) = -log I 0 c i A g i (A k )(A -A k ) + f i (A k ), (5)
where

g i (A k ) is the gradient of f i (•) at A k : g i (A k ) = -c i r 1 (c i A k ), (6)
with r 1 = I1 I0 and I 1 the modified first order Bessel function of the first kind.

Accordingly, the tangent majorant function

H L (•, A k ) of L is given as: H L (A, A k ) = N 2 A 2 σ 2 + N i=1 g i (A k )(A -A k ) + f i (A k ), (7)
and its gradient ∇H L (•, A k ) is expressed as :

∇H L (A, A k ) = N A σ 2 - N i=1 c i r 1 (c i A k ). (8)
The value of A k+1 canceling this gradient is directly derived as:

A k+1 = σ 2 N N i=1 c i r 1 (c i A k ). (9)
It can be seen that this update scheme corresponds to a classical gradient descent strategy:

A k+1 = A k - σ 2 N ∇L(A k ), (10)
with a fixed stepsize

α k = σ 2 N .
In the rest of this paper, this algorithm is referred to as MMχ as the Rician distribution is a special case of the non-central χ distribution.

IV. MM DESCENT ON A TRANSFORMED CRITERION

According to Proposition 3 in [START_REF] Idier | Properties of Fisher information for Rician distributions and consequences in MRI[END_REF], by introducing a change of variable y = A 2 , the transformed objective function becomes convex and unimodal with respect to the variable y. Therefore, we propose to construct a quadratic majorant function for the convex criterion and we will apply the descent scheme MM for the optimization with respect to y and deduce the value of A.

As illustrated by Figure 1, such scheme consists on constructing a non quadratic majorant in the initial space defined with respect to A which turns out to be more efficient than the quadratic surrogation around the concave zone. Indeed, after the transformation y = A 2 , the objective function is rewritten as:

L(y) = N y 2σ 2 - N i=1 log I 0 (c i √ y), (11)
Its gradient and Hessian are given by:

∇ L(y) = N 2σ 2 - N i=1 c i 2 √ y r 1 (c i √ y), (12)
∇ 2 L(y) = N i=1 c 2 i 4y r 2 1 (c i √ y) + 2 c i √ y r 1 (c i √ y) -1 . ( 13
)
Finally, as noted in section (2.2) of [START_REF] Laforgia | Some inequalities for modified Bessel functions[END_REF],

r 2 1 (x) + 2 x r 1 (x) -1 > 0 ∀x 0, (14)
and hence the positivity of the Hessian ∇ 2 L(y) can be deduced.

A. Construction of the majorant function

The proposed quadratic tangent majorant is given by:

H L (y, y k ) = 1 2 B k (y -y k ) 2 + ∇ L(y k )(y -y k ) + L(y k ), (15)
with a curvature defined by:

B k = 2 ∇ L(y k )y k -L(y k ) y 2 k . (16)
a) Discussion.: This curvature results from the specification to get a majorant function that passes through the same point as L at y = 0. This construction makes it possible to obtain a tangent quadratic majorant with the smallest curvature and thus allowing larger steps. The property of majorization can be deduced from the fact that the third derivative of the F is negative and therefore its curvature is decreasing. Demonstrations are omitted here due to the lack of space and we only show the proof of convexity of the majorant function H L (•, y k ).

Lemma 1: The function h(x) = xr 1 (x) -2 log I 0 (x) is negative for all x > 0.

Proof 1: Taking the derivative of this function

h ′ (x) = x -xr 1 (x) 2 -2r 1 (x)
and according to [START_REF] David | A tutorial on MM algorithms[END_REF], one can deduce that h(x) is a decreasing function. Therefore xr 1 (x) -2 log I 0 (x) < h(0) = 0, ∀x > 0, which completes the proof.

Theorem 1: The curvature B k of the tangent majorant function H(•, y k ) defined by Equation ( 16) is positive for any y k non-negative.

Proof 2: Note that B k (y) can be formulated as the first derivative of b(y) = 2 L(y) y , which can be expressed as

b(y) = N σ 2 -2 N i=1 log I 0 (c i √ y) y , with c i = Mi σ 2 . Its derivative is expressed by: b ′ (y) = - N i=1 c i √ y r 1 (c i √ y) -2 log I 0 (c i √ y) y 2
.

Lemma 1 allows to deduce the positivity of the curvature B k by setting x = c i √ y.

B. Majorant function minimization

The value of y k+1 canceling the gradient of H L (•, y k ) is expressed as:

y k+1 = y k -B k -1 ∇ L(y k ). ( 17
)
This algorithm is denoted as MM-χ 2 since it corresponds to applying the transformation y = A 2 . It can be noted that this update scheme is similar to a quasi-Newton algorithm where the pseudo-Hessian matrix is deduced from the Majorization principle.

V. PERFORMANCE ANALYSIS

This section is dedicated to an empirical analysis of the performances of the two MM optimization schemes in terms of convergence speed and the proposal of a new descent strategy combining the two algorithms.

A. Comparison of the two MM algorithms

In order to assess the performances for different noise levels, we first consider a simulation with the measured signal module to noise ratio (M/σ) varying from 0 to 8. The number of measurements is set to N = 1. A discussion on this choice will be conducted in section VI-A. The optimization is performed using the two MM algorithms. The measurement is taken as the initial value and the stopping condition is based on a gradient norm value less than 10 -5 .

Figure 2 shows the required number of iterations for the convergence of the two MM algorithms with different SNR values. We can observe that the necessary number of iterations for the convergence of MM-χ 2 is much smaller when M σ is less than a value in the order of 1.8, which demonstrates the interest of the variable change especially in the context of a low SNR. However, MM-χ algorithm requires fewer iterations when the SNR increases. This can be explained by the proximity of the Rician distribution to a Gaussian for high SNR.

B. Combination of the two MM algorithms

In order to obtain a more efficient algorithm for any value of the SNR, we firstly adopt a linear combination of these two descent algorithms either in the initial space A or in y. The optimization strategy will be referred to hereafter as a mixed descent whose main steps are summarized in Algorithm 1.

An optimized approach consists in selecting at each iteration the MM solution which induces the largest descent of the criterion (Algorithm 2). As illustrated by Figure 2, the proposed optimized descent and mixed descent make it possible to take advantage of the two algorithms by significantly reducing the number of iterations whatever the value of SNR. It is worth noting that the optimized descent turns out to take automatically the fewest number of iterations while the mixed descent results, for α = 0.5, in an intermediate convergence speed.

C. Performance analysis

A second experiment is performed to discuss the impact of the number of measurements on the convergence rate of the

Algorithm 1 Iterative mixed descent

Initialize A 0 = M Choose Descent in A or y Set α between (0, 1)

for k = 1, 2, 3, ... until convergence do Set y k = A 2 k
Compute g k and A k+1 according to Eq (6) and ( 9) Compute B k and y k+1 according to Eq ( 16) and ( 17)

switch Descent do case A Set D A = A k+1 -A k Set D y = √ y k+1 - √ y k A k+1 = A k + αD A + (1 -α)D y case y Set D A = A 2 k+1 -A 2 k Set D y = y k+1 -y k y k+1 = y k + α D A + (1 -α) D y A k+1 = √ y k+1 end for Algorithm 2 Iterative optimized descent Initialize A 0 = M for k = 1, 2, 3, ... until convergence do Set y k = A 2 k
Compute g k and A k+1 according to eq.( 6) and ( 9) Set L A = L(A k+1 ) as defined in eq.( 2) Compute B k and y k+1 according to eq.( 16) and ( 17)

Set L y = L( √ y k+1 ) if L y < L A then A k+1 = √ y k+1 end if end for
algorithms. A Monte-Carlo simulation with 10 5 realizations is performed with the true signal to noise ratio A/σ varies from 0 to 8 in the case of measurements number N = 1 or N = 10. Figure 3 illustrates the average number of iterations using the different variants of the MM algorithms. It can be seen that for both single and multiple measurements cases, the three proposed algorithms outperform the classical MM algorithm when the SNR is lower than a certain level and the Optimized descent improves significantly the convergence efficiency. Also, a rise in the number of iterations for A is noted for MM-χ when the number of measurements increases.

VI. APPLICATION TO MRI

This section is dedicated to the application to MRI data and a discussion on the relevance of the proposed algorithm.

A. Multiple or single measurements

Considering the dataset used for ML estimation in practical applications, two different raw images can be obtained according to the acquisition set-up. A common one is to collect N independent complex-valued images then convert them into N magnitude images. An alternative processing consists in accumulating the N complex-valued k-space images into a single observation by computing their average then processing the magnitude image [START_REF] Constantinides | Signalto-noise measurements in magnitude images from NMR phased arrays[END_REF]. The question is therefore whether the ML estimation should be performed by taking multiple measurements with single accumulation or a single measurement with multiple accumulations.

We compare these two approaches by an experiment with real MRI data of a cylinder phantom filled with oil. We collected ten repeated measurements of independent complexvalued observations, from which we can compute a number of magnitude measurements or alternatively taking the averaged magnitude image.

Figure 4 presents the mean absolute deviation and normalized mean square error of ML estimation given a number of independent measurements (from 1 to 10) compared with a given single image obtained by averaging the same number of accumulations. It can be observed that the second approach has a better performance in terms of deviation for some different number of accumulations. This demonstrates the superiority of ML estimation using the averaged image.

B. Application to a synthetic data sets

The proposed approach was evaluated using the freely available brain database [START_REF] Cocosco | Brainweb: Online interface to a 3D MRI simulated brain database[END_REF]. A set of realistic MRI data volumes with T 1 pulse sequence and 1mm slice thickness produced by the MRI simulator are taken. The noise variance is computed using a method of local noise variance distribution [START_REF] Aja-Fernández | Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach[END_REF]. In Figure 5, three different methods of ML estimation are applied to the noisy data: (b) represents the ML pixel-wise estimation where N = 1. As the ML estimation is originally designed to take multiple images and according to the results of section VI-A, a single averaged image is selected. We also perform two additional tests: (c) and (d) corresponds to ML estimation using 3 × 3 and 5 × 5 local neighborhoods where the number of measurements is 9 and 25 respectively. As mentioned by Santiago [START_REF] Aja-Fernández | Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach[END_REF], taking into account the local statistics for the estimation yields a better noise reduction but causes edge-blurring. Since this paper discusses the efficiency of different MM optimization strategies, we focus on the computational performance and not on the relevance of the ML estimation.

The MM-χ algorithm and the proposed three algorithms are evaluated with the same dataset and using the same stopping criterion in order to ensure getting identical solutions after convergence. To compare their convergence speed, their computation time (T), as well as the number of iterations (Iter) are reported in Tables I and II. It can be seen that the three proposed approaches substantially outperform MM-χ algorithm in terms of computational efficiency. Moreover, as the number of measurements increases, while MM-χ requires significantly more iterations, the Mixed descent remains a relatively moderate convergence speed while Optimized descent requires always the fewest iterations. In the case of pixelwise estimation where N = 1, we obtain the same result as in the previous simulation experiment, and for multiple measurements, the Mixed descent turns to be the most efficient in terms of the computation time. It is worth noting that even though for each iteration, the Mixed descent and Optimized descent are somehow more time-consuming as two criteria evaluations are involved, the combination strategy still results in less overall computation time due to fewer iterations.

VII. CONCLUSION

In this paper, we have investigated the problem of signal estimation for MRI with Rician noise. The main proposal is a novel MM framework based on a convex criterion which improves significantly the convergence efficiency in the context of a low SNR. The validation of its corresponding quadratic tangent majorant function has been demonstrated mathematically and empirically. Furthermore, we introduced two algorithms combining the two MM strategies in order to obtain the optimal performance.

From the simulation part and the real images analysis results, we can conclude that both algorithms have shown substantial improvement over MM-χ in terms of computational efficiency. A further theoretical analysis including the convergence rates of the algorithms is being investigated. The extension to penalized ML estimation is being explored in order to enhance the denoising performances.

Fig. 1 .

 1 Fig. 1. Illustration of the initial non-convex criterion and the transformed convex criterion and the tangent majorant functions in the initial space for (a) A k = 1.5, (b) A k = 7 and (b, d) after transformation y = A 2 .

Fig. 2 .

 2 Fig.2. The number of iterations of the MM descent algorithms for different levels of the signal to noise ratio in the case of a single measurement.

Fig. 3 .

 3 Fig. 3. Number of iterations of the MM descent algorithms for different original signal amplitude and (a) N = 1 measurement or (b) N = 10 measurements.

Fig. 4 .

 4 Fig. 4. The mean absolute deviation and normalized mean square error of ML estimation using a single measurement and multiple measurements.

Fig. 5 .

 5 Fig. 5. (a) Original noisy data. (b) ML pixel-wise estimation. (c) ML estimation using 3 × 3 neighborhoods (d) ML estimation using 5 × 5 neighborhoods

TABLE I THE

 I NECESSARY NUMBER OF ITERATIONS FOR THE FOUR MM ALGORITHMS FOR ML ESTIMATION FOR A SINGLE IMAGE USINGPIXEL-WISE, 3 × 3 AND 5 × 5 LOCAL NEIGHBORHOODS.

	Algorithm	Pixel-wise	3 × 3	5 × 5
	MM-χ	58	482	459
	MM-χ 2	47	47	46
	Optimised descent	9	11	11
	Mixed descent	14	14	14

TABLE II THE

 II COMPUTATION TIME (IN SECONDS) OF THE FOUR MM ALGORITHMS FOR ML ESTIMATION FOR SINGLE IMAGE USING PIXEL-WISE, 3 × 3 AND 5 × 5 LOCAL NEIGHBORHOODS.

		Pixel-wise	3 × 3	5 × 5
	MM-χ	5.8	308.1	839.7
	MM-χ 2	3.1	19.2	51.2
	Optimized descent	1.6	13.9	41.1
	Mixed descent	1.7	11.4	31.9

ACKNOWLEDGMENT

This work was funded by the GDR 720 ISIS in the framework of the 2016 call for projects.


Johan Lefeuvre

Saïd Moussaoui

Laurent Grosset

Anna-Luiza Mendes-Siqueira

Franck Delayens

Supervised Mixture Analysis and Source Detection from Multimodal Measurements

Keywords: Multimodal supervided spectral unmixing, Sensor fusion, Chemical mixture analysis

The aim of this paper is to present a method for source detection within unknown chemical mixtures using several measurement modalities. Contrary to the well studied case of single source detection, this approach enables simultaneous detection of multiple chemical components by exploiting the mixing coefficients resulting from supervised linear unmixing and thresholded non-negative least squares. The first contribution of this work is to propose an automated procedure to compute an optimized binary classifier rule for each component independently using a database of known mixtures. The second contribution is to propose a global decision rule based on the fusion of the multimodal decisions using weighting schemes such as those used in multiple classifier systems (MCS). A real database of Ion Mobiliy Mass Spectrometry (IMMS) data is used to evaluate the detection performance. An increase of the detection accuracy is reached using the multiple thresholds within the independent classifiers approach as compared to single modality detection.

Introduction

Source detection from spectral data is at the core of several applications of signal processing methods in physical sensing problems, such as chemical substance analysis [START_REF] Duarte | Source Separation in Chemical Analysis : Recent Achievements and Perspectives[END_REF], hyperspectral imaging [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF]. The concept of a source signal is defined as the spectral signature associated to a chemical component.

Classical approaches for source detection and identification are based on the recognition of some discriminant patterns or features of the sought sources [START_REF] Scharf | Matched subspace detectors[END_REF]. These features can be determined empirically [START_REF] Kopczynski | Peak modeling for ion mobility spectrometry measurements[END_REF][START_REF] Marczyk | Improving peak detection by gaussian mixture modeling of mass spectral signal[END_REF] or learned from a large-scale database [START_REF] Horsch | A detailed comparison of analysis processes for mcc-ims data in disease classification-automated methods can replace manual peak annotations[END_REF]. Often, the measured signal is interpreted as the linear combination of an unknown set of several sources. In such cases, multivariate analysis techniques [START_REF] Pomareda | Hard modeling multivariate curve resolution using lasso: application to ion mobility spectra[END_REF] such as independent component analysis and non-negative matrix factorization can be used. But, even if these blind separation methods yield estimates of the source and their relative abundances (mixing coefficients), their identification as a physically meaningful components is not guaranteed in all situations.

Alternatively, one can consider a known set of sought sources and estimate the mixing coefficients by a linear regression under non-negativity constraint. These coefficients are then used to make a decision on the presence or the absence of each source in the mixture. In this supervised linear unmixing framework, various methods can be applied, such as thresholded non negative least squares [START_REF] Slawski | Sparse recovery by thresholded non-negative least squares[END_REF], nonnegative orthogonal greedy algorithms [START_REF] Nguyen | Non-negative orthogonal greedy algorithms[END_REF] or constrained sparse regression [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]. The main challenge in these approaches is to estimate the sparsity level corresponding to the appropriate number of mixture components. For non-negative least-squares (NNLS), hard thresholding based methods may require a tuning of the threshold level [START_REF] Slawski | Sparse recovery by thresholded non-negative least squares[END_REF] or a setting of an appropriate stopping strategy of the iterative process [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF]. For greedy algorithms, the sparsity level should be defined manually or estimated automatically by adopting an adequate stopping rule during the decomposition. In this paper, we explore more specifically the hard thresholding based approach for which we propose two enhancements. The first one is to adopt a source dependent threshold which allows more accurate detection. The second proposal is to use a training database of known mixtures to determine an optimized detection threshold for each sought source separately.

The abundance of each source in an observed mixture data depends on the sensitivity of the sensing modality to this source. Consequently, a lack of sensitivity to some sources in one modality leads to a poor detection of these sources using this modality and conversely an accurate detection will be achieved from another modality more sensitive to these sources. Exploiting multi-modality will therefore enhance detectability of all the sources present in a mixture [START_REF] Fauvel | Decision fusion for the classification of urban remote sensing images[END_REF][START_REF] Jeon | Decision fusion approach for multitemporal classification[END_REF]. In the case of multiple measurements of the mixture recorded in the same modality, it has been shown that joint analysis of the data can enhance performance [START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF][START_REF] Tropp | Algorithms for simultaneous sparse approximation. part i: Greedy pursuit[END_REF]. In the case where different measurement modalities are available, an appropriate fusion strategy should be defined. The second contribution of this paper is to adopt a decision fusion method based on a multiple classifier system (MCS) [START_REF] Ho | Decision combination in multiple classifier systems[END_REF]. Finally, the proposed approaches of component-dependent thresholding and multimodal decision fusion are tested on a challenging example of chemical mixture analysis using Ion Mobility Mass Spectrometry (IMMS) data [START_REF] Kanu | Ion mobility-mass spectrometry[END_REF][START_REF] Szymańska | Chemometrics for ion mobility spectrometry data: recent advances and future prospects[END_REF], where spectral responses of the mixtures and the sources are recorded using four complementary measurement modalities (ionization modes).

Problem Statement

Let us consider the case of a single measurement modality. The measurement vector of the mixture is noted y ∈ R M , where M represents the number of samples provided by the sensor. Measurement vectors associated to N sources are gathered in a matrix S = [s 1 , ...s N ] ∈ R M ×N . The measurement model is assumed to be a linear mixing:

y = N i=1 a i s i + e, (1)
where the additive noise term e corresponds to measurement errors. Computing the vector of mixing coefficients a = [a 1 , ..., a N ] t ∈ R N can be efficiently done by solving the following problem:

â = arg min a∈C y -Sa 2 2 , (2)
where C denotes the constraint set of the coefficients. For the considered application, measured data correspond to mass spectrometry and ion mobility spectra and are assumed to follow a linear mixing model. The leas-squares problem above ( 2) is therefore solved under the constraint of non-negativity using a nonnegative least squares algorithm (NNLS) [START_REF] Lawson | Solving least squares problems[END_REF] or an interior-point least squares (IPLS) [START_REF] Chouzenoux | Fast constrained least squares spectral unmixing using primal-dual interior-point optimization[END_REF]. The estimated mixing coefficients in each measurement modality are then used to retrieve the detection vector d = [d 1 , ..., d N ] t ∈ {0, 1} N from a where for i ∈ 1, ..., N, d i = 1 when component i is present in the mixture and d i = 0 otherwise. In the case where L independent measurement modalities are available, for each modality l ∈ 1, ..., L, the measurement vector of the mixture is noted y (l) ∈ R M l and the spectra of the N sources are gathered in a matrix

S (l) = [s (l) 1 , ...s (l) N ] ∈ R M l ×N
where M l is the length of the data vector in the l-th modality. The observation model is then expressed as:

y (l) = N i=1 a (l) i s (l) i (3)
with linearly independent abundance vectors a (l) i . The detection of the components which are present in the mixture should therefore be realized by accounting for the values of the mixing coefficients in all modalities using a decision fusion strategy.

Figure 1 shows the proposed detection pipeline. It consists in firstly solving a non-negative linear regression problem in each modality and then deducing the binary detection vectors (d (1) , ..., d (L) ). The proposed method to compute a fused presented in Section 4. Section 5 details an application of the proposed detection method in the case of chemical mixture analysis using IMMS spectrometry data.

Detection from single modality measurements

A first step of the proposed approach consists in estimating the mixing coefficients (a i , for i ∈ 1, ..., N ) by solving Problem [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF]. A NNLS algorithm [START_REF] Lawson | Solving least squares problems[END_REF] is used for this purpose. Dedicated optimisation algorithms [START_REF] Chouzenoux | Fast constrained least squares spectral unmixing using primal-dual interior-point optimization[END_REF] can be used in order to account for additional constraints (such as sum-to-one). The detection of each mixture components from these mixing coefficient values is addressed as a binary classification problem between two states: presence or absence of each source.

Detection approach

The decision rule related to the two states leads to a binary classifier h θi (a i ) : R → {0, 1} according to:

h i (a i ; θ i ) = 1 if a i ≥ θ i 0 otherwise (4)
Our proposal consists in adopting thresholds θ i which are specific to each source in the considered measurement modality. It allows to account for the different sensitivities of the measurement modalities on the sought sources. The detection vector is then given by d i = h i (a i ; θ i ) for i ∈ 1, ..., N . The values of the thresholds θ i are defined by optimization of the detection performances on a training database of known mixtures. The best detection threshold values are chosen in such a way to reach a balance between sensitivity and specificity of the binary classifier.

Specification of the detection thresholds

Let Y = [y 1 , ..., y K ] ∈ R M ×K denote a matrix of K mixtures with a subset of known sources among a set of N sources. The detection performances for each binary classifier corresponding to the i-th source are evaluated in terms of true positive rate and false positive rate

TPR i (θ i ) = TP i TP i + FN i and FPR i (θ i ) = FP i FP i + TN i , for i = 1, . . . ,
(TPR i (θ i ), FPR i (θ i )) -(1, 0) 2 2 . (5)
Note that the ideal point (1, 0) corresponds to the detection of the source (sensitivity) when it is present in the mixtures and without any false positive detection (specificity).

Illustration on a real mixture data set

A real database of IMMS data from 85 lubricant mixtures is considered. Each mixture contains between 5 and 7 components among 20 possible ones. More details on this data set are given in Section 5. Figure 2-(a) gives the distance between the receiver operating characteristic (ROC) curve for the detection of different sources. It can be noted that the optimal detection thresholds for the considered sources (C6, C15 and C17) are different, which suggests to use a source dependent threshold. Moreover, the global threshold seems to be appropriate for C15 but it is not optimal for the two other sources (C6 and C17). shows a comparison between the average ROC curve obtained when applying a source dependent threshold and the ROC curve obtained by using the same threshold for all the sources. Both the Area Under Curve (AUC) and the average performance at the optimal point are higher in the case of source dependent thresholds.

Detection strategy in the multimodal case

This section addresses the detection in the case where L > 1 independent measurement modalities are available. To illustrate the relevance of this strategy let us compare the ROC curves, shown in Figure 3, obtained with two different sources C1 and C3 in the four modalities offered by the IMMS spectrometer. The ROC curves of modality M2 (to the left) yields the best performance and M3 the worst performance for source C1. In contrast, for component C3, the best detection is obtained by modality M1 and the worst detection is obtained with modality M4. This example illustrates the complementarity of the different modalities and the need for a decision fusion strategy for an accurate detection of the entire set of sources.

Weighting schemes

The fusion of the independent binary classifiers (h

i , . . . , h

(L)
i ) in the L modalities is performed using a linear combination rule according to

g i a (1) i , . . . a (L) i = L l=1 ω (l) i h (l) i a (l) i ; θ (l) i , (6) where θ (l)
i and ω (l)

i correspond to the detection threshold and the decision weight associated to the detection of i-th component in the l-th measurement modality. A resulting fused decision corresponds to a Multiple Classifier System [START_REF] Ho | Decision combination in multiple classifier systems[END_REF] defined subsequently as:

d fus i = 1 if g i a (1) i , . . . a (L) i ≥ 0.5 0 otherwise (7)
Depending on the values assigned to the decision weights, one can distinguish mainly three different fused classifiers. a) Majority Vote (MV). It consists of retaining the decision taken by the absolute majority of classifiers [START_REF] Boland | Majority systems and the condorcet jury theorem[END_REF]. It is defined as below:

ω (l) i = 1 L . (8)
However, the MV classifier doesn't account for the performance of the binary classifiers in the different measurement modalities.

b) Weighted Majority Vote (WMV). Weighted Majority Vote [START_REF] Kuncheva | Combining pattern classifiers: methods and algorithms[END_REF], [START_REF] Moreno-Seco | Comparison of classifier fusion methods for classification in pattern recognition tasks[END_REF] is a decision fusion rule in the case of independent classifiers and the weights are defined as:

ω (l) i = log µ (l) i 1 -µ (l) i (9)
where µ

(l)
i is the considered measure of the i-th classifier performance (distance to ideal point, global accuracy, balanced accuracy, etc.) in the l-th modality. c) Dynamic Classifier Selection (DCS). DCS [START_REF] Giacinto | Adaptive selection of image classifiers[END_REF] is a simple and powerful multiple classifier fusion strategy consisting in selecting the most efficient classifier for each source and discarding the others. The chosen decision weights in [START_REF] Duarte | Source Separation in Chemical Analysis : Recent Achievements and Perspectives[END_REF] depend on the performance index of the binary classifiers in each modality. The values of ω (l) i are set according to:

ω (l) i = 1 if µ (l) i = max k∈[[1,L]] {µ (k) i }. 0 otherwise. ( 10
)

Performance index

The weighting schemes for decision fusion are based on classifier performance. The most commonly used performance measure is the detection accuracy (ACC), defined as:

µ acc i = TP i + TN i TP i + TN i + FP i + FN i . ( 11
)
This index is commonly used for evaluation of classifiers with balanced occurrences of presence/absence of the sources.

In the case of our application, the source component are more often absent than present in the mixtures which leads trivial classifiers with high rejection rate to yield good global accuracy scores. This phenomenon called "curse of accuracy" [START_REF] Kubat | Addressing the curse of imbalanced training sets: one-sided selection[END_REF] is avoided by choosing a more adapted performance index such as the Mathews Correlation Coefficient (MCC) [START_REF] Powers | Evaluation: From precision, recall and F-factor to ROC, informedness, markedness & correlation[END_REF], which is defined as

µ mcc i = TP i × TN i -FP i × FN i (TP i + FP i )(TP i + FN i )(TN i + FP i )(TN i + FN i ) . ( 12
)
This index is considered as one of the best binary classifier performance metrics since it is not affected by class imbalances in the training set. An ideal classifier leads to µ mcc i = 1, a random classifier gives µ mcc i = 0 whether or not the training set is balanced, meanwhile a classifier systematically predicting the exact opposite of the ground truth will hace an MCC value µ mcc i = -1.

Application to mixture analysis by IMMS spectrometry

In this section, detection performances of the presented MCS are compared to the performance obtained separately on each modality on a supervised chemical mixture analysis by an IMMS spectrometer.

Mixture synthesis

A database of 85 mixture made from 20 different chemical components has been designed. Each mixture contains 5, 6 or 7 components randomly chosen from a set of 20 classical components involved in lubricant formulation [START_REF] Ponthus | Evaluating the multiple benefits offered by ion mobilitymass spectrometry in oil and petroleum analysis[END_REF]. Samples from each mixture and each source components have been analysed twice with an IMMS spectrometer [START_REF] Kanu | Ion mobility-mass spectrometry[END_REF]. This spectrometer ionizes the analyte sample to create a swarn of ions. In a second step, bidimensional maps corresponding to distribution of drift times and time of flight of the ions through two separation chambers are recorded. The distribution of drift time through the first chamber is called the ion mobility spectrum and the distribution of time of flight in the second chamber leads to the mass spectrum. Mass spectra and ion mobility spectra are considered as the spectral signatures of the analyte sample. Two distinct ionization modes have been used, therefore 4 one-dimensional spectra are recorded for each sample. Those modalities are called M1, M2, M3 and M4, corresponding to positive ionisation mass spectra, negative ionisation mass spectra, positive ion mobility spectra and negative ion mobility spectra. The spectral signatures of one of the mixture are presented in Figure 4-(a). Figure 4-(b) shows the spectral signatures of one component of this mixture.

Mixture analysis

The IPLS algorithm [START_REF] Chouzenoux | Fast constrained least squares spectral unmixing using primal-dual interior-point optimization[END_REF] is used for the mixing coefficient estimation under nonnegativity constraint. Among the 85 observed mixtures, a set of randomly chosen 60 mixtures, with balanced occurrences of each component, are retained for training the algorithms in terms of optimal threshold values setting and decision weights calculation. The remaining set of 25 observed mixtures are retained for the performance evaluation. This procedure is repeated with 60 independent realizations to get statistically robust global performance indexes.

Unimodal detection and majoruty vote

Global performance of each optimized unimodal classifier and the majority vote of the 4 modalities on the database are presented in Table 1. It can be noted that modality M2 seems to present the best global performance in terms of true positive rate, false negative rate and overall accuracy. It can be noted that a naive global fusion scheme such as MV leads to an improvement of specificity (lower TPR) but at the expense of sensitivity (higher FPR). This results suggests adopting alternative fusion strategies

Decision fusion for multimodal detection

The results of the application of the two fusion strategies based on binary classifier performances are reported in Table 2. Two weighting strategies based either on ACC and MCC indices are considered. For all the presented metrics, one can notice that these weighted methods (WMV and DCS) outperform the best of unimodal classifiers. More specifically, a significant improvement in ACC and MCC indices is achieved. An accuracy value of 87.9 % for the MCC based WMV. It notably reduced the rate of false positives (5.1%) while being less sensitive than the modality M2. However, the MCC based DCS classifier is the most sensitive one (82.3 %) and presents a moderate rate of false alarms (13.2 %).

Discussion on component detection

Figure 5 shows the detection results of two groups of sources. In the first group, the sources are well detected by positive ionization mode while in the second group, the sources are better detected using negative ionization mode. For each group of components, unimodal detection results are compared to the multimodal detection. In both cases, the performance of the fused decision is equivalent to the best of all the unimodal classifiers. The proposed framework has therefore been able to fully exploit the complementarity of the modalities.

Conclusion

The concrete problem of identifying chemical components in an unknown mixture from multi-model spectrometric data, while being an instance of the wellknown source separation problem, poses many challenges among which the sourcedependent response sensitivity and need to derive robust decision fusion strategies to combine information provided by the different modalities. It has been shown on real data that a more accurate detection is achieved by proposing a component-specific threshold in each modality and adopting a decision fusion scheme exploiting the detection performances in each modality, measured in a training database. We also noted that DCS classifiers tend to be very sensitive while WMV are more specific. Future works will be directed at investigating methods based on greedy sparse recovery, proposing adequate rules for the joint decomposition of the observed data with multiple modalities. Another perspective proposal of decision rules based on machine learning approaches and that will not requite the linear mixing model hypothesis.

Fig. 1 .

 1 Fig. 1. Multimodal detection pipeline

Fig. 2 .Figure 2 -

 22 Fig. 2. Influence of the threshold values on the detection performance

Fig. 3 .

 3 Fig. 3. ROC curves for source detection in four modalities

Fig. 4 .

 4 Fig. 4. Spectral responses of a mixture and one source in four measurement modalities

Fig. 5 .

 5 Fig. 5. Detection performance using measurements with two ionization modes

Table 1 .

 1 Performances of unimodal detection

	Modalities	Fusion
	Score M1 M2 M3 M4 MV
	TPR (%) 74.2 77.7 59.8 63.8 62.5
	FPR (%) 21.5 18.4 26.3 29.2 8.2
	MCC 0.50 0.57 0.31 0.32 0.57
	ACC (%) 77.4 80.8 69.7 69.3 83.9

Table 2 .

 2 Performance of fusion strategies based on either ACC or MCC measures.

		ACC	MCC
	Score WMV DCS WMV DCS
	TPR (%) 77.3 81.5 68.3 82.3
	FPR (%) 11.8 12.5 5.1	13.2
	MCC	0.65 0.67 0.62 0.67
	ACC (%) 85.5 86.0 87.9 85.6