Davinson Casta Ño-Cano
email: davinson.castanocano@cea.fr

Mathieu Grossard
email: mathieu.grossard@cea.fr

Arnaud Hubert
email: arnaud.hubert@femto-st.fr

Multi-axis Force Sensing with Pre-stressed Resonant Composite Plates: An Alternative to Strain Gauge Force Sensors

Industrial robots embedding multi-axis force sensors at the robot/environment interface presents numerous advantages in terms of safety, dexterity and collaborative perspectives. The key-point of these developments remains the availability of cheap but sufficiently precise multi-axis force sensors. This paper proposes a model-based approach to design a new alternative to commonly used strain gauge sensors. The principle of the device relies on pre-stress resonant composite plates where feedback control and measurement are achieved with piezoelectric transducers. The main originality of this work is that the force to be measured may present multi-axis components. Based on pre-stress and piezoelectric theories, a complete electromechanical model is proposed. This one is used during the design of a resonating composite Mindlin plate, embedding piezoelectric patches. It is shown that the effects of in-plane and out-of-plane external forces can be considered as pre-stress components. These ones, at the root of buckling phenomena, alter the stiffness of the structure and shift the plate resonance frequencies. Then, by solving the eigenvalue problem of the pre-stress vibrating structure, we can find the relationship between the natural frequencies of the structure and the externally applied multi-axis force. The proof of concept of this sensor is achieved on a case study. Finally, numerical results from both, home-made and commercial, finite element software demonstrates the interest of our approach to design integrated and inexpensive multi-axis force sensors solutions.

I. INTRODUCTION

Future challenges for industrial robotics concern mainly interactive, collaborative and safe robotics. In these perspectives, the ability to sense the arm-environment interactions is a crucial need. This is even the key-point for robots to evolve from purely repetitive behaviours to some degree of autonomy in unstructured surroundings and safe collaboration with humans [START_REF] Lamy | Human force amplification with industrial robot : Study of dynamic limitations[END_REF], [START_REF] Kosuge | Control of robot directly maneuvered by operator[END_REF]. These needs for accurate force/torque sensors have become the object of numerous research approaches in the literature [START_REF] Stefanescu | Handbook of Force Transducers: Principles and Components[END_REF].

Principles and technologies exploited for the sensor, as well as its functional integration into the sensor body device, remain key factors to define the performance and modality of the design (e.g. optical-based force sensors [START_REF] Chapuis | Design of a simple mri/fmri compatible force/torque sensor[END_REF], sensing capacitance array [START_REF] Sinden | A planar capacitive force sensor with six degrees of freedom[END_REF], modularized assembly of 1-D force sensors [START_REF] Wang | Toward force-based signature verification: A pen-type sensor and preliminary validation[END_REF], structures based on parallel Stewart kinematics [START_REF] Zhao | Performance analysis and optimization of sizable 6-axis force sensor based on stewart platform[END_REF], etc.). Currently, the majority of available force sensors employ strain gauges mounted on more or less conventional mechanisms (e.g. shear-type rosette strain gauges fixed onto torsion bars [START_REF] Fleming | Automotive torque measurement: A summary of seven different methods[END_REF] or Hollow Hexaform structures [START_REF] Aghili | Design of a hollow hexaform torque sensor for robot joints[END_REF] for torque sensors, Maltese cross shape for force sensors [START_REF] Uchiyama | A systematic design procedure to minimize a performance index for robot force sensors[END_REF], etc.), which can be fairly stiff and robust. However, strain-gauge-based sensors present drift problems and their manufacturing is process intensive (severe and precise requirements on the working, assembly and calibration of the integrated structure are needed). An attractive alternative to strain-gauge technology for force measurement are resonant-based sensors. Indeed, a great advantage lies in their easy way to measure frequency shift, and their reduced manufacturing cost. Moreover, resonant force sensors present advantages over static sensors in terms of noise immunity, high sensibility and good stability.

Because induced static stresses due to external forces slightly change the stiffness properties of the forced-vibrating structure, a slight frequency variation occurs. This results from the non-linearity of the strain field according to the external mechanical loads. Piezoelectric transducers, with their advantageous reversible (inverse and direct) electromechanical coupling effect, are usually employed in the force measurement principle based on integrated resonant sensors. Major studies currently found in literature about resonantoperating structures with multifunctional PZT patches rely on axially loaded beam elements [START_REF] Cheshmehdoost | Design and performance characteristics of an integrated high-capacity DETF-based force sensor[END_REF], [START_REF] Fukuzawa | Monolithically fabricated double-ended tuning-fork-based force sensor[END_REF]. Nevertheless, this design limits the force sensing capabilities to only one direction for each beam. This choice for one-axis sensor design can be explained mostly by drastic complications encountered by the design of multi-axis piezoelectric resonators. To exhibit interesting modal sensitivities to forces along many preferential sensing directions, a more complex geometry must be envisaged than simple beams.

The idea proposed in this paper is to exploit different resonant modes of a single composite plate as a way to measure the components of an externally applied force. A review of the literature shows that no previous paper have investigated the vibrations of composite plates with piezoelectric transducers, taking into account the effects of in-plane and out-of-plane forces. Force sensing capabilities of plates, as a multi-axis extension of beam's force sensing capabilities, have not been yet exploited. The mechanical pre-stress theory proves that forces in all directions can be measured [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]. The simple shape of a plate is easily manufactured, all the more that piezoelectric actuators and sensors patches can simply be bonded to the upper and lower faces of the plate. This simple structure avoids the complex spatial arrangement of 1D resonating beams [START_REF] Wang | Flexible building blocks: Modularized design concept of a six-axis force/torque sensor[END_REF] or the difficulty to accurately sticking many strain gauges on different faces of squared cross-section beams [START_REF] Stefanescu | Handbook of Force Transducers: Principles and Components[END_REF]. The ef-fects of in-plane (membrane effect) and out-of-plane external forces generate a pre-stress of the structure. These ones, at the root of buckling phenomena, alter the global stiffness of the structure and shift the plate resonance frequencies. A complete electromechanical model is established in this paper to design pre-stress and resonating composite Mindlin plate embedding piezoelectric transducers. The model makes use of the Hamilton's principle and energetic methods. Then, by solving the eigenvalue problem of the pre-stress vibrating structure, we can find the relationship between the natural frequencies of the structure and the externally applied multiaxis force. Due to the complexity of the computation, the eigenvalues of the resulting multiphysics problem is solved using an home-made finite element (FE) approximation.

The paper is organized as follows. In the next section, the model of the dynamic pre-stressed multilayer structure with piezoelectric transducers is established using variational principle. It explains how pre-stress phenomena influence the overall structure stiffness. In the third section, a FE formulation is proposed for Mindlin plate. Governing matrix equations of the full system are exposed. In the last section, a design approach for resonant force sensing is proposed through a case study. Our procedure to determine the applied multi-axis loads is finally validated by comparison with a commercial FE software.

II. MODELLING OF ACTIVE STRUCTURES WITH PRE-STRESS

Resonant force sensing is based on a mechanical structure whose natural frequencies and shapes are responsive to an external force (named F in the paper). Such force is considered as parts of the boundary conditions of the structure. The main idea lies on the fact that frequencies and modal shapes of a structure are strongly dependant on their boundary conditions (force or displacement).

A. Description of the Structure

The sensor device is a composite structure made with a structural plate and a set of piezoelectric sensors and actuators. The simplified design of this system is displayed in Fig. 1. Different boundary conditions are imposed (Dirichlet/essential conditions on displacements or Neumann/natural conditions on forces). Boundary forces are measured using the resonance frequencies shift. It has to be noticed that two kinds of forces are considered:

• The excitation force Fext appears at high frequency. Fext will be created with piezoelectric elements to excite the resonances of the structure. The resonances can be obtained by application of electric potential to the piezoelectric patch electrodes. • The to-be-measured force F appears at a relatively low frequency. The relation between an applied tobe-measured force F and the frequency shift of the structure can be understood via pre-stress loads. Indeed, in interactive robotics application and on the contrary of other forces Fext , the frequencies of to-be-measured force F are far below the resonant frequencies of the Description of the force sensor composite structure (in-plane dimensions length a and width b, thickness h). The system is made from a structural plate with collocated piezoelectric transducers for actuation and sensing. A set of forces F = { F1 , F2 , F3 } to be measured in different directions is applied at the edges of the plate (parts of boundary conditions).

V ε i j = εi j + εi j σ i j = σi j + σi j εi j V V (t) σi j F ůi t = 0 x 1 x 2 x 3
εi j σi j Fext ũi Fig. 2. The continuous volume V in a free state, is submitted to a preload F, that generates strain ε and stress σ into the body. The resulting configuration V is therefore seen as the initial configuration (t = 0), from which, this body is excited to its resonance in response to external forces Fext .

plate. In such a case, they can be considered as quasistatic forces that only affect the mean strain of the plate. This frequency separation creates two groups of external forces, low-frequency forces that will shift the resonance and high-frequencies forces that will put the plate into resonance. The main characteristics of the plate can then be split into two terms, quantities only due to pre-stress (shown with ˚), and quantities due to others sources (shown with ˜) [START_REF] Geradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]. For example, the strain can be split as ε i j = εi j + εi j where εi j is the pre-stress quasi-static contribution due to F and εi j a strain variation (due to Fext ) around the pre-stress configuration (Fig. 2).

According to the Green-Lagrange definition [START_REF] Landau | Theory of elasticity[END_REF], the strain ε i j in a continuous media is expressed in terms of its displacement u i as :

ε ik = 1 2 ∂ u i ∂ x k + ∂ u k ∂ x i ε (1) ik + 1 2 ∂ u l ∂ x i ∂ u l ∂ x k ε (2) ik (1)
Einstein's summation convention on repeated indexes is used in this paper. Scientific literature [START_REF] Geradin | Mechanical vibrations: theory and application to structural dynamics[END_REF] shows that in pre-stress structures, the main new contributions arise from Angles ˚Over-script for pre-stress terms ˜Over-script for excitation terms quadratic terms ε (2) ik that are often neglected in linear non-prestress mechanical problems. Therefore, the two terms ε (1) ik and ε (2) ik will be kept in the following sections. A brief description of each symbol used in the paper is shown in Table I.

B. Internal Energy for Passive Pre-stress Structures

Elastic energy density U m of continuous media results from an integration of stress σ along strain ε:

U m = ε i j 0 σ i j dε i j (2)
If stress and strain are split according to Fig. 2 (e.g. σ i j = σi j + σi j and ε i j = ε i j + ε i j ) and linear elastic constitutive law is assumed (Hooke's law: σ i j = c i jkl ε kl ), U m becomes:

U m = 1 2 c i jkl εi j εkl Ůk +c i jkl εi j εkl + 1 2 c i jkl εi j εkl (3)
In this expression, Ůk is a constant term. Other terms can be expanded considering entire definition of the Green-Lagrange strain (1) :

U m = 1 2 c i jkl εi j εkl Ůk + 1 2 c i jkl ε(1) i j ε(1) kl U k + σi j ε(1) i j U n + σi j ε(2) i j U g + c i jkl ε(1) i j ε(2) kl + 1 2 c i jkl ε(2) i j ε(2) kl U nl (4)
The second term U k is the classic first-order elastic strain energy density for non-pre-stress structure. As explained in [START_REF] Geradin | Mechanical vibrations: theory and application to structural dynamics[END_REF] and due to an internal equilibrium (force balance), the third term U n disappears when we consider the total energy of the system. The forth term, a second-order strain contribution, is known as the energy density of pre-stress [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF] or the geometric energy density U g [START_REF] Geradin | Mechanical vibrations: theory and application to structural dynamics[END_REF] (geometric because it does not depend on the material of the structure). This term will be exploited in the principle of our force sensor because it is directly related to the to-be-measured force F (through the term σ ). Finally the last higher-order terms (3 rd and 4 th order) named U nl , are required in the case of large displacements (i.e. for out-of-plane motions when small strain and small displacements theory is no longer valid).

C. Internal Energy for Piezoelectric Structures

Such force sensor device is therefore a sandwich compound of different materials (a passive structural plate where piezoelectric transducers are integrated, Fig. 1). Piezoelectric transducers are needed to generate the excitation Fext that puts the structural system in resonance.

• Piezoelectric actuators are dynamically excited to produce shear forces at the upper face of the plate. • Piezoelectric sensors permit to measure structural response in terms of surface strain. Where piezoelectric materials are present, the previous internal energy density has to be completed with the electric terms [START_REF] Tiersten | Hamilton's principle for linear piezoelectric media[END_REF]:

U = 1 2 (ε i j σ i j -E i D i ) (5)
If linear piezoelectric constitutive equations are assumed (superscript E stands for constant electric field, and ε stands for constant strain) [START_REF] Smith | Smart material systems: model development[END_REF], [START_REF] Landau | Electrodynamics of Continuous Media[END_REF]:

σ i j = c E i jkl ε kl -e ki j E k D i = e ikl ε kl + ε i j E j (6)
the internal energy density of piezoelectric parts becomes:

U = 1 2 c E i jkl ε i j ε kl U m -e ki j E k ε i j U em - 1 2 ε i j E i E j U e (7)
In this expression, we can distinguish the pure mechanical energy density contribution U m , the electromechanical contribution U em and the pure electrical contribution U e . It can be noticed that the mechanical contribution U m in [START_REF] Zhao | Performance analysis and optimization of sizable 6-axis force sensor based on stewart platform[END_REF], can be expanded for pre-stress structures, as previously (4).

D. Dynamic Modelling for Pre-stress Piezoelectric Structure

The dynamics of the sensor device is available using energetic method and variational principle. In this paper, this is done using the Hamilton's principle of least action and Lagrangian energetic function. The Lagrangian density of any electromechanical system is:

L (u i , φ ) = T * -U + W (8)
where T * is the kinetic co-energy density and W the work density done by external sources. For piezoelectric continuous media, the generalized displacements are mechanical displacements u i and electric potentials φ . The Hamilton's principle of least action states that:

δ A = δ t 2 t 1 V L dV dt = 0 (9)
for any virtual generalized displacements δ u i or δ φ . The solution of the related eigenvalue problem gives its natural frequencies and shapes (eigenvalues and eigenvectors).

In [START_REF] Fleming | Automotive torque measurement: A summary of seven different methods[END_REF], the excitation energy density W from external sources are of different natures. They come from electrical stimulus (electric potential φ and q electric charge density) or mechanical forces density f ext i , (Fig. 1):

δ W = f ext i δ u i -q δ φ (10)
In ( 8), the kinetic co-energy density of a continuous media is defined as (with ui = du i /dt):

T * = ρ ui ui /2 (11)
This model takes into account the mechanics of the structure and the electromechanical coupling phenomena due to the piezoelectric elements. The internal energy density U in (8) for a system considering piezoelectric and pre-stress effects, can therefore be written as a result of (4-7):

U = Ůk + U k + U n + U g + U nl + U em + U e ( 12
)
This energetic description allows us to describe the dynamic equations of the full sensor system. As stressed before, it is important to notice that the internal energy density [START_REF] Fukuzawa | Monolithically fabricated double-ended tuning-fork-based force sensor[END_REF], especially the geometric energy density U g , is a function of the to-be-measured force F. Therefore F is linked to σ which is linked to U g and finally to the variational principle δ A = 0. The to-be-measured force F is then directly embedded in the Euler-Lagrange equations and their modal shapes and frequencies. The solutions of this eigenvalue problem gives therefore a way to estimate the applied force F and then to the use of the device as a multi-axis force sensor.

Nevertheless, Euler-Lagrange equations for continuous media have analytical solutions only in the case of very simple geometries, or 1-D approximations, such as beams. Because of this lack of analytical solutions for plates with arbitrary boundary conditions, the model of this pre-stress plate will be solved numerically using a Finite Elements Method in the next section.

III. FINITE ELEMENT MODELLING

A. Multilayer Plate Element Description

The plate is modelled using Mindlin's theory, with the following assumptions:

• u 1 and u 2 displacements are linear functions of x 3 (linear displacements across the plate thickness, Fig. 3a). • Stress σ 33 , normal to the mid-plane, is negligible (null plane stress condition through the thickness: σ 33 = 0) if normal components of external forces Fext 3 are null. • Surfaces, normal to the mid-plane before deformations remains straights after, but not necessary normal to this plane (θ 1 and θ 2 not necessary equal to

∂ u 3 ∂ x 1 and ∂ u 3 ∂ x 2
, Fig. 3b and Fig. 3c) . For numerical and Finite Elements computations, Voigt's vector representation is commonly used instead of tensor one (Table II).

B. Finite Element Formulation

For a system with piezoelectric elements and considering pre-stressed initial conditions, the Hamilton's principle (9) (with Lagrangian (8) including kinetic [START_REF] Cheshmehdoost | Design and performance characteristics of an integrated high-capacity DETF-based force sensor[END_REF], external [START_REF] Uchiyama | A systematic design procedure to minimize a performance index for robot force sensors[END_REF] and potential ( 12) terms) can be written in Voigt's vector form as [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF], [START_REF] Piefort | Finite element modeling of piezoelectric structures[END_REF]:

0 = - V ρ {δ u} T { ü} -{δ ε (1) } T c E {ε (1) } -{δ ε (1) } T c E {ε (2) } -{δ ε (2) } T c E {ε (1) } -{δ ε (2) } T c E {ε (2) } -{δ ε (2) g } T [ σ ] {ε (2) g } + {δ ε} T [e] T {E} + {δ E} T [e] {ε} + {δ E} T [ ] {E} + {δ u} T F ext V dV + Sm {δ u} T F ext S dS - Se δ φ σ dS + {δ u} T F ext P -δ φ Q (13)
[ σ ] and {δ ε (2) g } are special matrix rearrangement terms described in [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF] to take advantage of Voigt's representation in the computation of U g . For clarity, superscript ˜is omitted in this section (superscript ˚is kept in the pre-stress matrix [ σ ]). Integrations are performed on pre-stressed volume V and surface S.

Finite element method is based on a finite-dimensional reduction with interpolation functions defined on nodes of a mesh of the continuum. In the case of an electromechanical system, we consider that displacement {u} and electric potential {φ } can be represented as a set of interpolation functions based on nodes values.

{u} = [N u ] {u i } (14a) {φ } = N φ {φ i } (14b)
The strain {ε} and electric {E} fields can be therefore computed using numerical differential operators [D] applied to the displacements and the electric potentials at the nodes:

{E} = D φ N φ {φ i } = B φ {φ i } (15a) {ε (2) g } = [D g ] [N u ] {u i } = [B g ] {u i } (15b) {ε (1) } = [D u ] [N u ] {u i } = [B u ] {u i } (15c) {ε (2) } = 5 ∑ ξ =1 {u i } T [B g ] T S ξ [B g ] {u i } = [B nl (u i )] {u i } (15d) x j x 3
x i (a) linear displacement across the thickness

x 1 x 2 x 3 θ 1 θ 2 (b) plate rotations u i u 3 θ i x j x 3
x i (c) x 3 ⊥ x i plane displacements Fig. 3. Representation of displacements of the plate according to the Mindlin's kinetic assumptions: three translations u 1 , u 2 and u 3 , and two rotations θ 1 and θ 2 . (with i, j = {1, 2} and i = j)

where, [D u ] is the numerical differential operator such that {ε (1) } = [D u ] {u} according to [START_REF] Lamy | Human force amplification with industrial robot : Study of dynamic limitations[END_REF]. D φ = -∇ is the numerical gradient operator, and [D g ] is a numerical differential operator described in [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]. The selection matrix S ξ rearranges each 2 nd order strain term of (1). As seen before, ( 13) including ( 14) and ( 15) must be valid for any virtual admissible variation of displacement δ u i or electric field δ φ i [START_REF] Piefort | Finite element modeling of piezoelectric structures[END_REF]. This leads to the matrix equations for each element of the system:

[M] { üi } + ([K g ] + [K uu ]) {u i } + K uφ {φ i } = f ext i K φ u {u i } + K φ φ {φ i } = {q i } (16)
where the mechanical stiffness [K uu ] contains a non-constant term as follows :

[K uu ] = [K l ] + [K nl ({u i })]
Matrix of mass, geometry (pre-stress), mechanical stiffness, electromechanical coupling and electric capacitance are respectively defined as:

[M] = V ρ [N u ] T [N u ] dV [K g ] = V [B g ] T [ σ ] [B g ] dV [K l ] = V [B u ] T c E [B u ] dV [K nl ] = V [B u ] T c E [B nl ] + [B nl ] T c E [B u ] + [B nl ] T c E [B nl ] dV K φ u T = K uφ = V [B u ] T [e] T B φ dV K φ φ = - V B φ T [ ε ] B φ dV (17)
External forces and electric charges are defined as:

f ext i = V [N u ] T F ext V dV + S f [N u ] T F ext S dS + [N u ] T F ext P {q i } = - Sq N φ T q dS -N φ T Q (18)
All the elements ( 16) can be assembled to obtain the matrix governing equations of the full system:

M 0 0 0 Ü Φ + K g + K UU K UΦ K ΦU K ΦΦ U Φ = F ext Q (19)
{F ext } is the vector of external forces applied to the structure and {Q} is the electric charges applied on the electrodes, the vector of unknowns is made by chaining the unknowns of each node as follows:

{U} = n ∑ i=1 {u i } and {Φ} = n ∑ i=1 {φ i } As [K g
] is a function of the measured force F, if a modal analysis of this system is performed, natural frequencies will be function of the to-be-measured force, according to our force sensor device principle.

Based on these analytical and numerical developments, and parts of codes from [START_REF] Ferreira | MATLAB codes for finite element analysis solids and structures[END_REF], this model has been implemented using the MATLAB software. A case study is proposed in the following section to see the relevance of the method.

IV. CASE STUDY OF A COMPOSITE PLATE

Simulations are conducted to find the structure resonant frequencies for multiple configurations of pre-loads. Relationships between force and frequencies are then built with a two-steps procedure, as detailed in Algorithm 1. Computation of the relationship between resonant frequencies and input force is achieved sequentially:

• The static analysis can be used to check the maximal stress generated by the to-be-measured force. Newton-Raphson method has been used to solve the nonlinearity due to K nl ( Ů) . • The modal analysis takes into account the strain field generated by F, and casts two kind of results that are used for different purposes: the eigenvalues (resonant frequencies) are used to set up the relationship between force and frequency. The eigenvectors (modal shape) are used to place optimally the sensors and actuators. This optimisation seeks to excite the structure with the minimum amount of energy and to get a high amplitude signal. In such a way, anti-nodes appear to be the best place to collocate the piezoelectric transducers.

Algorithm 1: Simulation steps to compute the relationship between resonant frequencies f i and input force F (steps are related with Fig. 2

([K g ] + [K l ] -λ [M]) ( V ) {U} = 0
with, λ = ω 2 and ω i = 2π f i end Output: Mapping F → f i : resonant frequencies as a function of the input force Due to the geometry of the structure, quadrangular shell elements (with 5 dofs) are used. The bonding between the passive structure and the piezoelectric elements is considered perfect (no glue layer is considered in this study). We assume that all piezoelectric patches are connected to electric ground (φ = 0), and the stiffness due to the electrical contribution will not be considered in the following.

In Fig. 4, we present a simple case of a composite plate made of a passive material in the middle part, covered with piezoelectric patches on both faces. Two adjacent edges are clamped whereas the to-be-measured force, possibly multidirectional, is applied at the opposite corner. The dimensions and mechanical properties of the sensing device are reported in Table III. The simulation results of our home-made FE MATLAB code are then compared with those of the commercial software COMSOL Multiphysics. As can be seen in Fig. 5, the first three resonance frequencies are clearly sensitive to the to-be-measured forces F and differently according to the three x 1 , x 2 and x 3 axis. It is therefore possible to use this device as a multi-directional force sensor.

Fig. 5(a) shows that the first resonant mode can be used to sense the external force along the three directions ( F1 , F2 , F3 ) whereas the second (respectively the third) is not adapted to track the external force along the x 2 direction (respectively x 1 direction): the resonant frequency shift is less sensitive to F1 (respectively F2 ) applied force.

Because of the symmetry of the structure in the x 3 direction, the orientation (+ or -) of the to-be-measured force F3 along the x 3 axis can not be distinguished by the shift frequency (only its amplitude is available). Indeed, the force in this direction will be seen as a tensile force irrespective of its orientation and this will increase the resonant frequency. But, this could certainly be overcome with an intelligent extraction algorithm or a fusion of additional rough static informations. Compared with COMSOL Multiphysics simulations, our home-made finite elements code gives similar results, especially for the x 1 and x 2 force axis. It can therefore be used with confidence for the design and calibration of such kind of force sensor devices.

V. CONCLUSIONS

This paper presents the concept of a new multi-axis force sensor. This one is based on pre-stress resonant composite plates and is an interesting and viable alternative to the commonly used strain gauge sensor. The case study proves that such kind of sensors are able to measure the three components of an externally applied force.

Moreover, as electric charge are constantly generated in the dynamic steady state working mode of this device, inexpensive electronic instrumentation can be used, on the contrary of strain gauge or static piezoelectric force sensor devices.

The next steps of this work will be the experimental validation of these results using a prototype currently under construction.

  Fig. 1.Description of the force sensor composite structure (in-plane dimensions length a and width b, thickness h). The system is made from a structural plate with collocated piezoelectric transducers for actuation and sensing. A set of forces F = { F1 , F2 , F3 } to be measured in different directions is applied at the edges of the plate (parts of boundary conditions).

Fig. 4 .

 4 Fig.[START_REF] Chapuis | Design of a simple mri/fmri compatible force/torque sensor[END_REF]. Case study: composite rectangular plate with clamped conditions in two adjacent edges and a point of force F applied at the free corner.

  First resonance around 834Hz, mode (1,1). Second resonance around 2487Hz, mode (1,2).

Fig. 5 .

 5 Fig. 5. Relationship between the three force components and the first three resonant frequencies for the case study structure (includes the forcefree modal shape).

TABLE I NOMENCLATURE

 I

	Symbol	Description
	c	Elasticity modulus
	D	Electrical displacement field
	E	Electric field
	e	Piezoelectric material matrix
	F V , F S , F P	Body, Surface, Punctual forces
	K	Stiffness matrix
	K g	Geometric matrix
	M	Mass matrix
	Q	Electric charge
	q	Electric charge density
	U	Internal energy density
	V, V	Body free of any load, and pre-stressed body
		Permittivity
	ε	Mechanical strain
	σ	Mechanical stress
	ρ	Mass density
	φ	Electric potential
	θ x , θ y

TABLE II NOMENCLATURE

 II FOR Voigt'S VECTOR QUANTITIES

		Mechanical		Electrical
	{F}	Force	q	Electric Charge
	{u}	Displacement	φ	Electric Potential
	{σ }	Stress	{D} Electric Displacement
	{ε}	Strain	{E}	Electric Field

  )

	Input: Boundary conditions.
	for F do
	1. Static Analysis : Pre-stressing the structure
	Input: Force vector F
	• Calculation of the displacement field Ů as
	solution of : Ů = [K l ] + K nl Ů -1 (V )	F
	2. Modal Analysis : Measuring the force via a
	frequency shift
	Input: Pre-stressed displacement field Ů
	• Calculation of the geometric matrix [K g ( σ )]
	• strain ε calculation (15) from Ů
	• stress σ calculation via Hooke's law
	• Geometric stiffness [K g ] (17)
	• Computation of resonant frequencies f i :


Naresh Marturi
email: naresh.marturi@femto-st.fr

Sounkalo Dembélé

Nadine Piat

Scanning Electron Microscope Image Signal-to-Noise Ratio Monitoring for Micro-Nanomanipulation

Keywords: Scanning electron microscope, signal-to-noise ratio, image quality estimation, performance evaluation, micronanomanipulation

As an imaging system, scanning electron microscope (SEM) performs an important role in autonomous micro-nanomanipulation applications. When it comes to the sub micrometer range and at high scanning speeds, the images produced by the SEM are noisy and need to be evaluated or corrected beforehand. In this article, the quality of images produced by a tungsten gun SEM has been evaluated by quantifying the level of image signal-to-noise ratio (SNR). In order to determine the SNR, an efficient and online monitoring method is developed based on the nonlinear filtering using a single image. Using this method, the quality of images produced by a tungsten gun SEM is monitored at different experimental conditions. The derived results demonstrate the developed method's efficiency in SNR quantification and illustrate the imaging quality evolution in SEM.

Introduction

Vision-based control, also known as visual servoing, refers to the use of image information in order to control the motion of a robotic device. This technique uses the images acquired from a vision sensor to extract the information reflecting relative spatial position between the camera and the robot [START_REF] Hutchinson | A tutorial on visual servo control[END_REF][START_REF] Chaumette | Visual servo control. I. Basic approaches[END_REF]. In the recent years, it has become very popular especially in the field of robotic manipulation due to its ability in producing accurate and robust control. This has encouraged many researchers to extend this technique to micro-nanoscale in order to solve various problems associated with micro-nanomanipulation. Basically at this small scale (< 10 µm), manipulations are often performed using a scanning electron microscope (SEM), a transmission electron microscope (TEM) or an atomic force microscope (AFM). With its ability of producing images at high magnification in real-time, a SEM is always favored as an imaging tool for automatic micro-nanomanipulation applications.

In general, visual servoing control strategies are velocity-based and require a minimum frame rate of about 20 Hz. This corresponds to a very high scan rate compared to the normal operation of SEM. Moreover, to perform a vision-based micromanipulation task using SEM, one main requirement is that the captured images are to be noiseless and the details present in the images are high enough to be exploitative. Unfortunately, the noise level is critical in a SEM, especially at high acquisition rates. Each step during the image formation process adds noise which arise the difficulty in understanding and handling the final image noise. Practically, noise variation can be determined by acquiring images with varying scan time per pixel. If this time increases, the level of signal increases as the beam will spend more time at the particular location. However, with an increase in the scan time, the time to capture one frame increases which is not an acceptable constraint for high speed visual servoing applications. So, in order to perform a visual servoing-based micromanipulation task, it is mainly required to monitor the acquired final image quality until the overall task is accomplished. This monitoring process has been given the major attention in this paper.

In order to estimate the image quality, the primary requirement is to identify the type of final image noise. [START_REF] Reimer | Scanning Electron Microscopy: Physics of Image Formation and Microanalysis[END_REF] provided a detailed description of SEM and pointed out the fact that the noise added due to the primary electron emission is shot noise and grow in proportion to the useful signal. It is added due to the fluctuations in the count of emitted primary beam electrons. This type of noise follows Poisson distribution. For SEMs with thermionic guns, this shot noise is more (Sim et al., 2004b). [START_REF] Timischl | A statistical model of signal-noise in scanning electron microscopy[END_REF] have described five sources of noise in a SEM: primary emission, secondary emission, scintillator, photocathode, and photomultiplier. They have showed every noise correspond to quantum fluctuations and can be modeled by Poisson statistics. However, many works assume that the noise added by the detection system is negligible (Sim et al., 2004b;[START_REF] Goldstein | Scanning electron microscopy and X-ray microanalysis[END_REF]. [START_REF] Mulapudi | Is SEM noise Gaussian[END_REF] have showed empirically that in a thermionic gun SEM, final image noise follows Gaussian statistics but not Poisson statistics. In order to generate virtual SEM images, [START_REF] Cizmar | Simulated SEM images for resolution measurement[END_REF] have considered that the final image noise is an addition of a Poisson distribution representing primary emission and a Gaussian distribution representing the others types of noise in the SEM (secondary emission and electronics). Sim et al., (2004a;[START_REF] Sim | Image noise cross-correlation for signal-to-noise ratio estimation in scanning electron microscope images[END_REF]2012) have studied the final noise in SEM images. They assumed this noise is additive white noise and have developed many methods for estimating image signal-to-noise ratio (SNR) based on image cross-correlation. Since these previous investigations create ambiguity in selecting the type of noise, in this work, we first study the final image noise.

To assess the quality of acquired frames, image SNR has been selected as an evaluation parameter in this work. SNR is a commonly used term in the field of signal processing to differentiate or to estimate the level of noise present in a recorded signal [START_REF] Bose | Digital signal and image processing[END_REF]. In general terms, it is defined as the ratio of desired signal over background noise and is given by Equation (1).

10

(1) So far in the literature, many research works have proposed various methods using two SEM micrographs for the computation of SNR using auto-correlation technique [START_REF] Frank | Signal-to-noise ratio of electron micrographs obtained by cross correlation[END_REF][START_REF] Erasmus | Reduction of noise in TV rate electron microscope images by digital filtering[END_REF]. The main drawback associated with these methods is that the acquired images are to be perfectly aligned and are obtained by scanning the same specimen area. In addition, they require long processing time and are complex to integrate and use with real time systems. This work has been performed to monitor the image quality during an automatic vision-based micronanomanipulation task. In this paper, we present a simple and fast SNR estimation method using a single image. An estimate of the image noise to use with the developed method has been computed using non-linear spatial filtering. As the overall processing time is less, the developed method can be used with real time applications. Using the developed method, the image SNR of a tungsten filament gun SEM has been monitored during various experimental conditions that includes varying scan times, magnifications, operating time and device focus as well as during an automatic nanopositioning task.

Materials and Methods

Experimental set-up

The main set-up used for the experiments is shown in Figure 1. It consists of a JEOL JSM 820 SEM equipped with tungsten filament gun, a conventional Everhart-Thornley SE detector, an image acquisition system (DISS5 from point electronic GmbH), a SEM control computer (Intel Pentium 4, CPU 2.24 GHz and 512 MB of RAM) and a work computer (Intel Core 2 Duo, CPU 3.16 GHz, and 3.25 GB of RAM). The accelerating voltage for the SEM varies from 0.3 kV to 30 kV and the magnification varies from 10x to 100,000x. The maximum allowable electronic working distance is 50 mm. The communication between the two computers is accomplished using a client-server model. The server program runs from the SEM control computer that is interfaced with SEM control electronics and the image acquisition system. The image acquisition system is mainly responsible for sending the scan parameters to SEM and to acquire the data coming from secondary electron detector. Later this data is amplified, digitized by a 4 x 12 bit A/D converter (DISS5) and saved as an image in the server computer. The acquired image is then transferred to the client over TCP/IP where the process of monitoring takes place. The SEM chamber is also equipped with three degrees of freedom (xyz) piezoelectric positioning platform (TRITOR 100 from Piezosystem Jena GmbH). It can provide a motion up to 100 µm on all the three axes with a resolution of 0.2 nm in open loop. Different axes of this piezo positioning platform are controlled by the NV 40/3 piezo voltage amplifier. For visual servoing tasks, the platform control voltage values are transferred to the NV 40/3 from the work computer via RS-232 (serial port).

Study of final image noise

As mentioned, noise sources are numerous in SEM imaging and each contribute independently towards the final image formation. At the first hand, there are random statistical fluctuations in the primary and secondary electron emissions which produce shot noises independent from one pixel to another and depends on the material being used for imaging. Generally this type of noise satisfies the Poisson statistics:

! "! $ % (2)
where, λ the number of occurrences of the events and µ is the mean signal level. However, if the mean number of events is large enough (µ >10), the Gaussian distribution is a good approximation for the Poisson distribution.

On the other hand, we have also noise added by the detection and recording system. Apart from this, noise is also added by the photomultipliers where they emit some signal even in the absence of original signal. This noise is due to the thermal fluctuations in the photocathode. Upon considering all these noises, the final image noise n(x, y) at the pixel (x, y) is a result of superposition of all these partial noises. Under the central limit theorem this final noise can be approximated with Gaussian statistics [START_REF] Bovik | Handbook of image and video processing[END_REF].

Developed method for SNR quantification

From the above study, we assume that the final acquired image &, ( of silicon microparts sample shown in Figure 2

As an initial step, the contrast of the image is adjusted using histogram equalization. This is optional as most of the present day microscopes come with an integrated module to perform this task during image acquisition.

Next, the median filter mask of appropriate size is selected based on the amount of blurring and the total amount of processing time taken. The amount of blurring has been computed by calculating the mean squared error (MSE) between the original image and the filtered image * &, ( given by Equation ( 4) [START_REF] Bose | Digital signal and image processing[END_REF].

From the noise image, it is clear that the slight edges of microparts (can be treated as signal) are also appeared in the noise image. Since, estimating the noise using a single image is under constrained, it is difficult to estimate the accurate noise level by filtering. However, many works have pointed out this problem and proposed some solutions to tackle this issue [START_REF] Tai | A fast method for image noise estimation using laplacian operator and adaptive edge detection[END_REF]. Since this work is mainly aimed at performing quick monitoring of image quality to assist the final visual servoing task, we do not consider estimating the accurate amount of noise.

Finally, using the estimated noise free image * &, ( and noise image &, ( , the SNR is computed using Equation ( 7) obtained from Equation (1).

>

? @ A B,C ? D B,C E (7)
where, F √H8I78 J$, F * K,L and F K,L are the standard deviations of filtered and noise images, respectively.

Experimental samples and conditions

All the experimental codes are programmed in C++ and are executed from the work computer. The micro-nanostructures used for various experiments are fabricated in the micro-manufacturing center at FEMTO-ST.

The sample used for the noise study, scan time, magnification and focus experiments is a gold on silicon microscale calibration pattern [START_REF] Marturi | Vision and visual servoing for nanomanipulation and nanocharacterization using scanning electron microscope[END_REF] containing multiple chess grid regions with squares ranging from 25 µm down to 1 µm per side. All the experiments are conducted using an aperture size of 50 µm, 10 kV accelerating voltage, at 1000x magnification. The magnification test is a completely automated process where the magnification values are transferred to the DISS5 device from work computer via image server. The magnification is digitally changed by varying the size of the scan length on the sample. The drift observed at high magnifications has been corrected automatically using the method presented by Marturi et al. (2013b). For focus experiment, the DISS5 image acquisition system used for this work provides a simple control for the focus by linking the objective lens focal length with a series of focus steps (i.e. each focus step modifies the focal length to get a focused image). The focus steps are varied automatically and more details can be found in (Marturi et al., 2013a). The images are acquired with a scan time of 360 ns for both magnification and focus experiments.

For monitoring the image quality in real-time the sample used is a silicon wafer containing the microstructures of dimensions 6 µm x 6 µm x 2.5 µm (Figure 2). They were designed to be assembled by their respective U-shaped grooves (without using any external adhesives) to form three dimensional structures [START_REF] Dembélé | Gluing free assembly of an advanced 3D structure using visual servoing[END_REF]. These type of 3D structures can be used for building complex micro-electro-mechanical systems (MEMS). For this test, the magnification used is 5000x, accelerating voltage is 10 kV, scan time of 360 ns and aperture size of 30 µm. The sample used with automatic positioning experiment is the silicon micropart (Figure 3) of dimensions 10 µm x 100 µm x 20 µm. This experiment has been conducted in order to characterize the surface stiffness of the micropart. For this test the images are acquired using an accelerating voltage of 10kV at 1000x magnification. For all the experiments, the SEM secondary electron images of dimensions 512 x 512 pixels are used.

Results

Study of the final image noise

This test has been conducted to study the final image noise with fast scan rates. It is performed by acquiring the images of calibration pattern at different scan times and by computing the noise distribution within a single region (Figure 4(A) and 4(B)). The scan times used for this test are: 180 ns (least possible with our system) and 360 ns. Two plain image regions are selected for each image where one region contains only gold and the other contain only silicon. Figure 4(C) to 4(F) show the intensity histograms along with the fitted distributions (Gaussian and Poisson) for gold and silicon regions at different scan times.

Filter size selection

The artificial noise free and noisy images (Figures 5(A) and 5(B)) of gold on carbon generated by Artimagen library [START_REF] Cizmar | Simulated SEM images for resolution measurement[END_REF]) are considered for this purpose. It is a library provided by the National Institute of Standards and Technology (NIST) to produce artificial SEM images with varying noise. The additive Gaussian noise with a standard deviation equal to 0.07 has been added to generate a noisy image. The MSE between original image and noisy image is found to be 110.30. The computed MSE values along with the total processing time taken using different filter sizes are summarized in Table I. The plots shown in the Figure 6(A) to 6(D) are the intensity values of the original image and filtered image along the horizontal line shown in the Figure 5(A).

Precision testing of the developed approach

The precision of the developed method is tested by using a noise free image (Figure 5(A)) and by corrupting it with a known amount of noise (SNR). Later the SNR values of the corrupted images are computed using the proposed method and are compared with the known values. Table II summarizes the computed values along with the percent error between the filtered image and the original image.

Quality monitoring with respect to scan time

Initial tests are performed to evaluate the image quality with a change in the scan time. It is a well-known fact that the SNR increases with increase in scan time. This experiment has been performed in order to test and support this statement. Nine different images of gold on carbon sample (at 30,000x magnification) and calibration pattern (at 1000x magnification) are acquired with varying scan times. Figure 7 shows the result of monitoring process.

Quality monitoring with respect to magnification

The second tests are performed to monitor the SNR evolution with change in magnification. For this test, images of the calibration pattern are acquired for the magnifications ranging from 250x to 1000x with a step of 10. For each magnification one image is acquired and the SNR is computed directly at the time of acquisition.

Figure 8 shows the evolution of the SNR with respect to magnification.

Quality monitoring with respect to focus

In this test, the images of the calibration pattern are acquired for a series of focus steps (change in focus) and the image SNR has been monitored. Figure 9 shows the variation of SNR along with image sharpness with change in the focus steps. The maximum of the sharpness curve provides the best focus position. For this test the image sharpness has been computed using normalized variance sharpness function given by Equation ( 8) (Marturi et al., 2013a).

(

)

2 1 1 ( , ) M N S g x y MN µ µ = - ∑∑ (8)
where, M and N are the image dimensions, µ is the mean intensity of the image ( , ) g x y given by Equation ( 9). 1 ( , )

M N g x y MN µ = ∑∑ (9)

Quality monitoring with respect to time

In this test, the SNR quantification and image quality monitoring is performed with instrument operational time. The total operational time selected is 45 minutes. As the overall operating time is larger, the specimen surface is connected to the mounting plate using a copper conductive tape to remove the accumulated charge on the sample surface. Total time taken to acquire, transfer and process a single image is about 400 milliseconds. The obtained SNR values during this test are plotted and shown in the Figure 10. The output data has been fitted using polynomial fitting in order to see the actual response.

Quality monitoring with respect to time during nanopositioning

The characterization of an optoelectronic microstructure is now considered. The structure is arranged horizontally on the platform, the vertically disposed probe comes into contact with the structure to estimate its stiffness. The structure must be positioned more accurately and quickly as possible under the probe. The best solution is to use the images provided by the SEM to achieve the control of the platform. But if the signal level is not high enough, so the SNR high enough, it will not be possible to measure the information necessary for accomplishing the task. Monitoring the SNR of images should allow avoiding using the images that do not have the level of SNR required. This required level will be determined after performing visual feature extraction tests at different experimental conditions and is not discussed in this paper. The experiment starts after stabilization about 40 minutes after the ignition of the SEM. Figure 11(A) and 11(B) show the reference and initial images. Figure 11(C) to 11(E) shows the error during the task and Figure 11(F) shows the final error or finish of the positioning task. The positioning error is computed and compensated by the visual servoing control strategy [START_REF] Marturi | Vision and visual servoing for nanomanipulation and nanocharacterization using scanning electron microscope[END_REF]. The corresponding evolution of the SNR of the images for each iteration is shown in Figure 12.

Discussion

After analyzing the intensity histograms and fitted distributions (Figure 4), it is clear that any of the two distribution curves (Poisson and Gaussian) do not fit correctly with the real pixel intensity distributions. However, Gaussian distribution seems to provide better fit than Poisson distribution which can be clearly seen with the intensity histograms especially with silicon regions.

From Table I, it can be seen that filter size of 5x5 establishes a good compromise between accuracy and time of SNR computation with the current setup. With another experimental setup, the choice may vary. Upon testing the precision of the developed method in estimating the original SNR, the mean variation of the error is comparatively less for 5x5 filter (Table 2) and the percent error is less than 10 in most of the cases (with a filter size of 5 x 5). However, with high amount of noise in the image (especially for SNR < 12 dB) the percent error is more than 10. This is mainly due to the inability of median filtering in the presence of high amount of noise. In fact, a median filter is able to remove noise in the images only if the noise pixels are less than one half of its neighborhood. This is the reason why a 3x3 filter provides poor error and a 7 x 7 filter provides better error than a 5 x 5 filter (for some conditions). However, from the performed experiments with our system, it has been found that the SNR level is more than 14 dB in all cases.

Monitoring the SNR with respect to scan time has quantified a result known qualitatively, namely that increasing the dwell time leads to an increase in the SNR of images. The obtained results make this clear i.e., the image quality increases with increase in the dwell time. This phenomenon can be explained as that the increase in dwell time reduces the overall raster rate and increases the amount of time the beam spends at a particular location. Because of this, the electron count recorded at the detector increases subsequently improving the image quality. Besides, to control a robot using visual feedback high frequency acquisition rate is required that is possible only by using a lowest dwell time. We must therefore ensure that the dwell time selected provides a sufficient level of SNR.

The evolution of the SNR with respect to magnification is a nonuniform decreasing function. With an increase in the magnification, since the size of the monitor (image) is fixed, the size of the scan area is narrowed and the surface corresponding to scan a pixel reduces. Therefore, the area of the gold region that actively participate in generating the dominant quantity of original signal reduces due to which the SNR reduces. This phenomena mostly depends on the material of the sample used for imaging.

The results obtained by monitoring the imaging quality with respect to focus show that the SNR level varies simultaneously with image sharpness. This phenomenon can be explained in two ways. The primary one is, with a change in the focus step (working distance), the size of beam interacting with the sample surface varies. At the point of in-focus, the current density is high resulting in more signal from the sample improving the SNR. The second one is, with a change in the focus, the level of original signal details present in the image varies. At the point out-of-focus, noise level is predominantly more than that of the signal level, e.g. at the initial steps of the sharpness curve, the images are completely dark containing no image details other than noise.

There are two aspects in the evolution of image SNR over time: a random aspect and a dynamic aspect. The first aspect is due to the inherently random nature of the noise as well as electronic signal. The secondary electron emission from the interaction of the primary electrons with the electrons of the atomic cloud on a sample surface is a random phenomenon. The second aspect is characterized by a rapid increase of SNR to a peak for a time of about 500 s, followed by a rapid decay to a minimum for a time of about 750 s, and an irregular growth up to stabilization for a time of the order of 2500 s. This time corresponds to the time taken by the filament to reach a specific temperature to produce stable number of electrons. As a result it is recommended to wait at least 40 minutes from the ignition of JSM 820 before starting an experiment. The successful implementation of nanopositioning using image feedback (Figure 11) confirms the validity of this recommendation.

Conclusion

The SEM is by far an important tool for producing static images with high resolution. Today, the trend is to use the SEM as a vision sensor for performing robotic manipulation of microstructures. In this case the observed scene is dynamic and requires fast image acquisition. But electronic imaging is by far known for its sensitivity to the scanning frequency i.e., a high frequency of acquisition leads to low quality images. In this case, it is necessary to analyze the image quality beforehand. In this paper, we have developed a method to quantify the level of image SNR. It is based on image non-linear filtering in the spatial domain and can be used in real time for quality monitoring. Using the developed method, image quality of a tungsten gun SEM has been monitored with respect to the following parameters: scan time, magnification, focus and operational time. The main lesson that can be drawn from these experiments is the instability of the SNR: it varies randomly and greatly depending on time as well as magnification. The analysis of the average change over time revealed two phases. In the first phase, the ratio increases and decreases suddenly. It corresponds to the rise in temperature of the filament. The second phase is a slow growing of the ratio that corresponds to a temperature stability of the filament over time. The great merit of this work done allows the quantification of the instability of SEM imaging. In the future, the method will be applied to two other types of SEMs (a LaB6 gun SEM and a FEG SEM) to study the evolution of the SNR of their images.

LIST OF FIGURES

LIST OF TABLES

  (A) is formed by the addition of additive white noise &, ( to the original noise free image &, ( .

  and N are the image dimensions. In order to accelerate the overall process, a fast sorting algorithm is used in developing the median filter. Now, the goal is to filter &, ( by applying the median filter to find a closest estimate * &, ( of the noise free image &, ( . Thus * &, ( 6$ 78 9 :, ; < (5) where, :, ; are the local filter neighborhood pixel coordinates selected from the acquired image. Figure 2(B) shows the filtered image * &, ( of Figure 2(A). In turn this filtered image is subtracted from the original image to produce a noise image &, ( given by Equation (6). Inverted version of the noise image (black regions correspond to noise) is shown in Figure 2(C). &, ( &, ( 1 * &, (

Figure 1 :

 1 Figure 1: Experimental set-up.

Figure 2 :

 2 Figure 2: Images of silicon microparts demonstrating the developed method. (A) Final acquired image (B) filtered image using a filter size of 5x5 (C) noise image.

Figure 3 :

 3 Figure 3: Silicon microparts used for positioning experiments.

Figure 4 :

 4 Figure 4: Images of gold on silicon calibration pattern acquired at (A) 180 ns/pixel (B) 360 ns/pixel. Intensity histograms with approximated distributions for (C) and (E) gold regions (D) and (F) silicon regions for both images respectively.

Figure 5 :

 5 Figure 5: Artificially generated (A) noise free and (B) noisy images of gold on carbon using Artimagen library.

Figure 6 :

 6 Figure 6: Plots of intensities along the horizontal line shown in Figure 5(A) after filtering using the filter size of (A) 3x3 (B) 5x5 (C) 7x7 (D) 9x9.

Figure 7 :

 7 Figure 7: Evolution of SNR with respect to scan time.

Figure 8 :

 8 Figure 8: Evolution of SNR with respect to magnification.

Figure 9 :

 9 Figure 9: SNR and image sharpness variation with respect to focus steps.

Figure 10 :

 10 Figure 10: Evolution of SNR with the device operational time.

Figure 11 :

 11 Figure 11: Series of images depicting the nanopositioning task. (A) Desired position. (B) Initial image. (C) -(E) Errors at different positions. (F) Final error.

Figure 12 :

 12 Figure 12: Evolution of SNR during nanopositioning task. Error positions shown in Figure 11 are pointed with the boxes.

  FIGURES

  Figure 7

Figure 12

 12 Figure 12

Table I :

 I Computed MSE values and processing time with various filter sizes.

Table II :

 II Known and computed SNR values using the developed method for artificial gold on carbon image with different filter sizes.

This work is conducted with a financial support from the project NANOROBUST (ANR-11-NANO-006) funded by the Agence Nationale de la Recherche (ANR), France. It is also performed in the framework of the Labex ACTION (ANR-11-LABX-01-01) and the Equipex ROBOTEX (ANR-10-EQPX-44-01) projects.


Mohamed Taha Chikhaoui

Kanty Rabenorosoa

Nicolas Andreff

Kinematic Modeling

Keywords: Continuum robot, Active cannula, Electro-Active Polymers (EAP), Kinematic modeling

HAL is

Introduction

Continuum robots are still enthralling researchers' interests, almost half a century after the first early prototype: the "Tensor Arm" of [START_REF] Anderson | Tensor arm manipulator design[END_REF]. Exceptional usefulness of continuum robots appears in applications where it is restrictive to have joints and stiff links. They have the potential to suffer localized damage while still maintaining a healthy degree of functionality [START_REF] Gravagne | Kinematic transformations for remotely-actuated planar continuum robots[END_REF]. Moreover, continuum robots are able to navigate through complex anatomy. They may be considered as part of robots for MIS (Minimally Invasive Surgery) and NOTES (Natural Orifice Transluminal Endoscopic Surgery). Recent medical continuum robots include laparoscopic application, laser manipulators, catheters and micro-endoscopes, summarized in [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots : A review[END_REF]. Concentric tubes developed by [START_REF] Sears | A steerable needle technology using curved concentric tubes[END_REF][START_REF] Webster | Toward active cannulas : Miniature snake-like surgical robots[END_REF] for endonasal skull base surgery are a major contribution to continuum robots category. Such manipulators are constructed from three precurved telescopic concentric tubes that have small diameters (less than 3mm). They are actuated at their base by translation and axial rotation of each tube, and the overlapping builds the final shape. Concentric tubes are the starting point of our study. We aim to change tubes curvatures by the means of embedded actuators, and thereby, provide additional degrees of freedom to the system. Monitoring the curvature and bending each tube in different directions are indeed expected to enhance the manipulator performances. Our considered manipulator is presented in Figure 1a, showing a CAD design of a curved tube that holds a laser tool at its end-effector. The bending is performed with Electro-Active Polymer (EAP) electrodes coated around the tube based on [START_REF] Shoa | Conducting polymer actuator driven catheter -overview and applications[END_REF] work and sketched in Figure 1b. We present in the next section the kinematic modeling of an active cannula robot, starting from a standard approach to achieve modeling of an EAP actuated continuum manipulator. In section 3, both models will be used to analyze a part of the reachable workspace. Then we can establish a theoretical comparison between the existing configurations and our design, using also performance indices to confirm the manipulability improvement. Finally, in section 4, we will conclude and present several future challenges that still need to be achieved in this field.

Kinematic modeling

Before we can start modeling differential kinematics, a few assumptions need to be set up in the standard approach. This concerns the description used in previous work [START_REF] Webster | Toward active cannulas : Miniature snake-like surgical robots[END_REF]. Furthermore, the modeling of a continuum robot with EAP-based embedded actuators will be described.

Standard approach

Arc parameters

The standard approach [START_REF] Webster | Toward active cannulas : Miniature snake-like surgical robots[END_REF] is based on a configuration of n concentric tubes presented in Figure 2a. Each tube i ∈ {1..n} is made of a straight part (S i ) and a precurved distal one (C i ) and can translate (by ρ i ) and rotate (by θ i ) with respect to the z-axis of the robot base. Thus, depending on the translation of each tube, the concentric assembly can be decomposed into m successive links, defined by the concentric overlaps of straight parts, precurved parts or nothing (e.g. for three tubes, C 3 /S 2 /S 1 or C 3 /C 2 / / 0). Each link is modeled by an arc of a circle (constant curvature assumption), described by three parameters: its length ℓ j , its curvature κ j which is the inverse of the radius of curvature r j , and the angle φ j of the so-called equilibrium plane containing the arc (Figure 2b).

Depending on the overlapping of the n tubes, the curvature of link j ∈ {1, ..., m} is given by:

κ j = κ 2 x j + κ 2 y j with κ x j = ∑ n i=1 E i I i κ i, j cosθ i, j ∑ n i=1 E i I i , κ y j = ∑ n i=1 E i I i κ i, j sinθ i, j ∑ n i=1 E i I i (1)
where E i is the elastic modulus, I i is the cross sectional moment of inertia, κ i, j is the intrinsic curvature of the i th tube in the j th link and θ i, j denotes the i th tube angle about the j th link frame zaxis. Finally, the equilibrium plane angle is given by:

φ j = arctan(κ y j /κ x j ) (2)
a) b)

Fig. 2: Active cannula description: a) three concentric tube configuration [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots : A review[END_REF], and b) arc parameter description [START_REF] Webster | Mechanics of precurved-tube continuum robots[END_REF].

Specific mapping and independent mapping

In [START_REF] Webster | Mechanics of precurved-tube continuum robots[END_REF], three spaces are defined:

• the actuator space: {ρ i , θ i |i ∈ {1..n}} • the configuration space: κ j , φ j , ℓ j | j ∈ {1..m} • the task space: SE(3)

Two space transformations are thus defined:

1. The specific mapping from the actuator space to the configuration space (actuator dynamics). This mapping totally depends on the actuation of the tubes. 2. The independent mapping from the configuration space to the task space (forward kinematics). This mapping is the same for all concentric tube architectures, satisfying the assumption of constant curvature links and can be generically modeled.

Forward kinematics can be accomplished in a variety of ways: through Denavit-Hartenberg parameters [START_REF] Hannan | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF], Frenet-Serret frames [START_REF] Hannan | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF], integral formulation [START_REF] Chirikjian | A modal approach to hyper-redundant manipulator kinematics[END_REF], and exponential coordinates [START_REF] Sears | A steerable needle technology using curved concentric tubes[END_REF] and [START_REF] Webster | Mechanics of precurved-tube continuum robots[END_REF]. Using the latter convention, the transformation T j from link j -1 to link j decomposes into a rotation of center r j = [1/κ j , 0, 0] T about the y axis by α j and a rotation about the z axis by φ j :

T j = R z (φ j ) 0 0 1 R y (α j ) p j 0 1 (3)
where α j = κ j ℓ j and p j = [r j (1cosα j ), 0, r j sinα j ] T .

Differential kinematics

To compute constant curvature kinematics of a multi-section tube, one must compute single section tube kinematics. For the brevity of this paper, the computation is omitted but details can be found in [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots : A review[END_REF]. The velocity of link j with respect to link j -1 is given by:

V j =         cos∆ φ j (cos(κ j ℓ j ) -1)/κ 2 j 0 0 sin∆ φ j (cos(κ j ℓ j ) -1)/κ 2 j 0 0 -(sin(κ j ℓ j ) -κ j ℓ j )/κ 2 j 0 1 -ℓ j sin∆ φ j 0 -κ j sin∆ φ j ℓ j cos∆ φ j 0 κ j cos∆ φ j 0 1 0         J j   κ j ∆ φ j l j   where ∆ φ j = φ j -φ j-1 (4)
Using the adjoint transformation introduced by [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF], the full independent kinematic Jacobian can be deduced from the individual ones:

J indep = J 0 Ad T 0 J 1 Ad T 01 J 2 . . . Ad T 0(m-1) J m (5)
where T 0 j =T 0 T 1 . . . T j is the j th transformation matrix at the j th link and

Ad T = R R[t] × 0 R (6)
with R and t the rotation and translation component of T and [t] × the skewsymmetric matrix associated to the vector cross-product by t. Thus, determining the number of links in a configuration is a preliminary task to understand full Jacobian matrix dimension. For a configuration with three totally curved concentric tubes, we obtain a 6 by 9 matrix, as there are 3 links. The specific kinematic Jacobian maps actuator derivatives ρi θi T into arc parameters derivatives κ j ∆ φ j l j T . In the three totally curved concentric tube configuration, the full specific kinematic Jacobian J spec is a 3 by 2 matrix. Consequently, this configuration is a non-holonomic robot: the whole space can be reached but, in a given state, only a subset of velocity directions is achievable.

Modeling of EAP actuated concentric tubes

Adding an embedded actuation to the previous configuration is beneficial as it provides a direct control of the intrinsic curvatures κ i, j of each tube, whereas the standard approach only takes into account constant intrinsic precurvatures. This adds one component per tube in the actuator space (ρ i , θ i , κ i ) but does not change the two other spaces. Therefore, the independent mapping and the independent kinematic Jacobian are the same as those described above, and only the specific mapping and specific kinematic Jacobian need to be derived. The curvature of a tube made of Electro-Active Polymer (Figure 1) follows a linear law in terms of voltage according to the relation explained in [START_REF] Shoa | Conducting polymer actuator driven catheter -overview and applications[END_REF]:

κ i, j = 1/r i, j = C PPY V i , (7)
where C PPY is considered as the Polypyrrole actuation constant and V i is the applied voltage.

By appropriate low-level control, we foresee to be able to servo κ i, j , the intrinsic precurvature of tube i in link j, to a desired value. Thereby, the specific mapping will have the same expression as in the standard approach, with only a change in the inputs (precurvatures). To compute the dependent Jacobian, for three totally curved concentric tubes, we first need to determine the number of links. With such configuration, there are three links: three tubes for the first link, two tubes for the second link, and one tube for the third link as shown in Figure 3a.

Differentiating (1) with respect to κi, j and θi, j yields:

κ j = 1 κ 2 x j + κ 2 y j κ x j κ y j A j B j C j D j κin, j θin, j (8)
where κin,j = [ κ1, j κ2, j κ3, j ] T , θin,j = [ θ1, j θ2, j θ3, j ] T , κi = C PPY . Vi and

A j = E 1 I 1 cosθ 1, j Σ i E 1 I 1 E 2 I 2 cosθ 2, j Σ i E 2 I 2 E 3 I 3 cosθ 3, j Σ i E 3 I 3 (9) B j = - E 1 I 1 κ 1, j sinθ 1, j Σ i E 1 I 1 - E 2 I 2 κ 2, j sinθ 2, j Σ i E 2 I 2 - E 3 I 3 κ 3, j sinθ 3, j Σ i E 2 I 2 (10) C j = E 1 I 1 sinθ 1, j Σ i E 1 I 1 E 2 I 2 sinθ 2, j Σ i E 2 I 2 E 3 I 3 sinθ 3, j Σ i E 3 I 3 ( 11
)
D j = E 1 I 1 κ 1, j cosθ 1, j Σ i E 1 I 1 E 2 I 2 κ 2, j cosθ 2, j Σ i E 2 I 2 E 3 I 3 κ 3, j cosθ 2, j Σ i E 2 I 2 (12)
Similarly, differentiating (2) yields:

φ j = 1 κ 2 x j + κ 2 y j -κ y j κ x j A j B j C j D j κin, j θin, j (13)
For the j th link, the corresponding arc length derivative corresponds to l j = ρ j . This leads to the specific kinematic Jacobian at the j th link:

  κ j φ j l j   = 1 κ 2 x j + κ 2 y j    κ x j A j + κ y j C j κ x j B j + κ y j D j 0 -κ y j A j + κ x j C j -κ y j B j + κ x j D j 0 0 0 κ 2 x j + κ 2 y j      κin, j θin, j ρ j   ( 14
)
The above specific kinematic Jacobian is a square matrix and does not contain any structural non-holonomic constraint. Moreover, it is a 3 by 3 matrix. Consequently, only two tubes should be enough to reach any pose in the workspace. However, keeping three tubes provides us with a redundancy which is highly recommended, especially in medical and surgical manipulators.

3 Performance analysis

Workspace analysis

Workspace means the reachable zones of a manipulator end-effector. In traditional robotics, one can obtain the workspace by inverting the direct geometrical model. However, in continuum robotics, converting the effect of curvature change or bending angle into movements is significantly more complex. We restrict our analysis to planar movements for an intermediate step. This is sufficient to illustrate the key advantages of changing an additional arc parameter: the curvature. We take into account a simple configuration case: three tubes with the same insertion angle α i that are not rotated about their z axis. The total robot length is 45mm and the tube outer diameters are respectively 3.05mm, 1.45mm and 0.72mm. The only possibility for a precurved concentric tube configuration to reach additional zones in the xz plan is to combine tube translations. This would be controlled by the precurvatures already defined and thus, provides a reduced freedom to the end-effector. Moreover, path-controlling this movement would be noticeably challenging. However, monitoring all the curvatures yields more options to the manipulator to sweep even more space. Figure 3a shows that changing the second tube intrinsic curvature κ 2,2 in the second link, provides an additional reachable zone, without even changing first and third tube translation or insertion angle. The workspace was generated with a translation sampling of 3mm, and a second tube curvature change sampling of 20m -1 . The workspace previously reachable by the existing active cannula type robot is hold, and an additional set of movements is achieved each time one of the curvatures is changed. It is proven by both workspaces superimposition shown in Figure 3b.

Performance indices

This analysis is based on the full Jacobian J robot singular values σ i . Thus, we need a singular-value-decomposition (SVD) of the matrix. This study is based on the three most significant performance indices: manipulability, isotropy, and condition number. Mathematical exact definitions can be found in [START_REF] Angeles | Performance evaluation and design criteria[END_REF] and their expressions are :

Manip = ∏ i ρ i , Isotropy = ∏ i ρ i ∑ i ρ i , Condition number = ρ max ρ min (15)
As shown in Figure 4, monitoring the tube pre-curvature has a significant effect. Firstly, we notice that manipulability is enhanced when the tube bends beyond κ 2, j =50m -1 (Figure 4a). Otherwise, it decreases as the tube is straightened to a linear configuration. Nevertheless, the manipulator is able to reach additional zones, which was impossible with constant curvature tubes. The same phenomena is observed for the isotropy index. Straightening the tube draws the manipulator closer from a singular configuration as the isotropy reaches zero. In both cases, we notice that the isotropy measure is very low. It is due to the robot architecture that undeniably does not allow velocities in all directions similarly. Observing the third curve, the singular position is confirmed when κ 2, j is near zero, with a condition number close to infinity (Figure 4c). Beyond κ 2, j =50m -1 , the Jacobian matrix is well-conditioned. The straight position of the manipulator is a singular configuration; thus, it is more challenging to achieve velocities. The manipulability indices are increasing with the curvatures. It is owing to the easiness to generate radial movements in a bent position.

Conclusion

In this paper, a novel EAP-based actuation technique of an active micro-endoscope was briefly described. Demonstrating the benefits of monitoring the tubes curvatures in contrast with a constant curvature existing configuration was the main contribution of this paper. On the one hand, this has been proven in terms of dependent Jacobian analysis which changes from non-holonomic to holonomic. On the other hand, performances of both manipulators have been compared through a part of the workspace as well as via the three most significant performance index evaluation. It has proven that the variable curvature improves the continuum robot performances.

For the future works, additional mechanical constraints as shearing and torsion have to be included. Moreover, other tube materials that allow more flexibility need to be explored. Another challenge is to improve the actuator design: modifying the electrodes patterned along the tubes would enable more bending directions and would provide more degrees of freedom to the manipulator.

Fig. 1 :

 1 Fig. 1: Design of the considered continuum robot. a) CAD design of the tubes, b) Illustration of an EAP-based embedded actuation [8].

Fig. 3 :

 3 Fig. 3: End-effector workspace is marked by blue '+' for constant curvature active cannulas and by red 'o' for flexible continuum robot. a) Schematic description of tube translation and second tube curvature κ 2, j change from 50 to 10m -1 . b) Workspaces superimposition in the xz plane.

Fig. 4 :

 4 Fig. 4: Performance index variation according to second tube precurvature variation below and beyond κ 2, j =50m -1 fixed for the configuration in the standard approach in Section 2: a) Manipulability, b) Isotropy, and c) Condition number.

Acknowledgment

This work has been supported by the Labex ACTION project (contract ANR-11-LABX-01-01) and by µRALP, the EC FP7 ICT Collaborative Project no. 288663.