A Ferrand
email: adrien.ferrand@cea.fr

M Darmon

S Chatillon
email: sylvain.chatillon@cea.fr

M Deschamps
email: m.deschamps@i2m.u-bordeaux1.fr

Modeling of ray paths of head waves on irregular interfaces in TOFD inspection for NDE

Keywords: head wave, ray tracing algorithm, irregular surface, wave diffraction, non-destructive evaluation

The TOFD (Time Of Flight Diffraction) technique is a classical ultrasonic inspection method used in ultrasonic non-destructive evaluation (NDE). This inspection technique is based on an arrangement of two probes of opposite beam directions and allows a precise positioning and a quantitative evaluation of the size of cracks contained in the inspected material thanks to their edges diffraction echoes. Among the typical phenomena arising for such an arrangement, head waves, which propagate along the specimen surface and are chronologically the first waves reaching the receiver, are notably observed. Head wave propagation on planar surfaces in TOFD configurations is well known. However, realistic inspection configurations often involve components with irregular surfaces, like steel excavated specimens.

Surface irregularity is responsible for numerous effects on the scattering of bulk waves, causing the melting of surface and bulk mechanisms in the head wave propagation. In order to extend the classical ray approach on these complex cases, a generic algorithm of ray tracing between interface points (GIRT) has been designed. With respect to time of flight minimization (i.e. the generalized Fermat's principle), ray paths can be computed by GIRT for different natures of waves scattered by the complex surfaces or by flaws. The head wave fronts computed by GIRT are notably in good agreement with FEM simulated results. This algorithm, based on pure kinematic analysis of waves propagation, represents a first step in the future development of a complete ray theory for head waves simulation on irregular interfaces.

Introduction

Head waves are well known phenomena occurring in seismic inspections: they correspond to the first wave received by seismic probes, giving so their names. Typical seismic inspections deal with configurations involving two media separated by a planar interface. Specific wave theories have been developed, all considering head waves on planar surfaces as critically refracted waves. These modeling approaches are based on an asymptotic integral evaluation [START_REF] Brekhovskikh | Waves in layered media[END_REF], the asymptotic ray theory [START_REF] Friedrichs | Geometrical acoustics. II. Diffraction, reflection, and refraction of a weak spherical or cylindrical shock at a plane interface[END_REF][START_REF] Červený | Theory of seismic head waves[END_REF] or the generalized ray theory (Cagniard-De Hoop method) [START_REF] Cagniard | Reflection and Refraction of Progressive Seismic Waves[END_REF][START_REF] Hoop | A modification of Cagniard's method for solving seismic pulse problems[END_REF]. It is shown that the head wave on planar interfaces corresponds to the combination of both waves travelling along one interface side with the corresponding material sound velocity (compulsory more than that of the incident wave), and bulk waves radiated from the interface towards one material and propagating at the critical angle. The waves propagating along the interface are composed of the well-known lateral waves [START_REF] Brekhovskikh | Waves in layered media[END_REF][START_REF] Friedrichs | Geometrical acoustics. II. Diffraction, reflection, and refraction of a weak spherical or cylindrical shock at a plane interface[END_REF][START_REF] Červený | Theory of seismic head waves[END_REF][START_REF] Cagniard | Reflection and Refraction of Progressive Seismic Waves[END_REF][START_REF] Hoop | A modification of Cagniard's method for solving seismic pulse problems[END_REF][START_REF] Uberall | Surface waves in acoustics[END_REF] and of the Goodier-Bishop type waves [START_REF] Goodier | A note on critical reflections of elastic waves at free surfaces[END_REF] for the case of a finite beam impacting the interface.

Head waves on planar surfaces are also observed in NDE by using specific techniques and similar configurations than the seismic ones. Among the various NDE inspections techniques, the Time of Flight Diffraction (TOFD) technique allows a precise and quantitative evaluation of surface breaking or embedded cracks. This technique uses two probes moving on the inspected material surface with a fixed Probe Center Spacing (PCS). The TOFD configuration is particularly well adapted to generate and receive head wave signals [START_REF] Mahaut | TOFD inspection simulation using Civa software[END_REF][START_REF] Chaffaï | Simulation Tools for TOFD Inspection in CIVA Software[END_REF][START_REF]Assessment, Management of ageing of major nuclear power plant components important to safety primary piping in PWRs[END_REF] (sometimes used for surface breaking cracks detection).

The CIVA software platform for NDE simulation, developed at CEA/LIST, enables a complete simulation of TOFD inspections: indeed different semi-analytical models [START_REF] Mahaut | TOFD inspection simulation using Civa software[END_REF][START_REF] Darmon | Simulation of disoriented flaws in a TOFD technique configuration using GTD approach[END_REF][START_REF] Zernov | A refinement of the Kirchhoff approximation to the scattered elastic fields[END_REF] have been integrated to model all the echoes observed in a TOFD configuration (head wave, flaw edges diffractions, back wall echo…). The developed model for head waves simulation is based on the asymptotic ray theory [START_REF] Červený | Theory of seismic head waves[END_REF] which has also been used to improve the modeling of flaw corner echoes [START_REF] Huet | Modelling of corner echo ultrasonic inspection with bulk and creeping waves[END_REF][START_REF] Mahaut | Modeling of corner echo in UT inspection combining bulk and head waves effect[END_REF]. Nevertheless, the existing head wave simulation can only be applied for a planar specimen composed of an isotropic medium.

TOFD and other pitch-catch configurations are also classically used on non-planar interfaces both in NDE (inspections of cylinders, nozzles [START_REF] Mahaut | TOFD inspection simulation using Civa software[END_REF]..) and in the seismic domain: indeed in works carried out by Zhou and Chen [START_REF] Zhou | Ray path of head waves with irregular interfaces[END_REF], head wave travel times have been calculated by numerical simulations on some irregular surfaces for some seismic inspections.

It has been shown that the head wave propagation is more complex on irregular surfaces than on planar surfaces. Indeed, whereas the head wave on a planar interface is due to the critical refraction of the incident wave and includes a surface wave component propagating along the planar surface (called lateral wave), the head wave propagation on irregular surfaces implies propagations in the material bulk. The received head wave signal on irregular surfaces has a different amplitude than in the planar case and a travel time compatible with the hypothesis of propagation in the material bulk. Analytical modelings of the head wave propagation around the Earth have also been proposed in several seismic studies. Indeed, using an integral representation of the solution, the received head wave signal has been evaluated for slightly curved irregularities [START_REF] Lerche | A mean-field solution of the reflection of a spherical acoustic wave from a rough interface[END_REF][START_REF] Lerche | On the reflection of acoustic waves from a slightly curved interface[END_REF] of surface parametric equation ( ) ² / 2 F x x D = with large values of the radius of curvature D compared to the surface size, and for spherically symmetric radially heterogeneous media [START_REF] Hill | Critically refracted waves in a spherically symmetric radially heterogeneous earth model[END_REF]. These studies confirms conclusions given in [START_REF] Zhou | Ray path of head waves with irregular interfaces[END_REF] concerning the time of flight and the amplitude of the head wave.

In the present study, complete interpretation and simulation of the head wave propagation mechanisms for irregular surfaces are proposed. This approach leads to the development of a generic algorithm of ray tracing between interface points (called in the following for simplification Generic Interface Ray Tracing/GIRT) which is able to compute the travelling path of head waves in specimens inspected by TOFD technique. Based on an adaptation to the NDE domain of seismic ray tracing algorithms, GIRT is a new method for solving the twopoints ray tracing problem, which is the calculation of the true ray between two defined points: the source and observation points.

Ray-tracing techniques are often employed for modeling wave propagation in electromagnetism, optics, seismology … Traditional ray tracing algorithms are the shooting method [START_REF] Julian | Three-dimensional seismic ray tracing[END_REF] and the bending method [START_REF] Wang | Tomographic inversion of reflection seismic amplitude data for velocity variation[END_REF][START_REF] Wang | Seismic amplitude inversion in reflection tomography[END_REF]. The shooting method, well adapted for field propagation modeling, fixes the source point, takes initial ray direction and then finds the location of the observation point; in the view of two-points rays tracing, an iterative procedure for finding the starting direction can be used to yield the desired observation point. Ray bending is a variational approach which perturbs some initial estimate of the ray path, the ends being fixed, until the true ray which minimizes the travel time is found. The two previous traditional methods are heavy in time computation for solving the two-points ray tracing problem and not efficiently adapted to deal with diffracted waves and shadow areas in which head waves around irregular surfaces propagate.

To improve efficiency of the shooting method on the diffracted waves calculation, seismic ray tracing algorithms based on a grid based scheme were developed more recently : wave front construction [START_REF] Vinje | Traveltime and amplitude estimation using wavefront construction[END_REF], the finite difference Eikonal equation solver [START_REF] Vidale | Finite-difference calculation of travel times[END_REF][START_REF] Vidale | Finite-difference calculation of traveltimes in three dimensions[END_REF][START_REF] Koketsu | Finite difference traveltime calculation for head waves travelling along an irregular interface[END_REF][START_REF] Rawlinson | Multiple reflection and transmission phases in complex layered media using a multistage fast marching method[END_REF][START_REF] Kool | A practical grid-based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media[END_REF] which also generally needs large computation time for tracing rays and Shortest Path Method (SPM) algorithms [START_REF] Nakanishi | A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure[END_REF][START_REF] Moser | Shortest path calculation of seismic rays[END_REF][START_REF] Cao | Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm[END_REF][START_REF] Van Avendonk | Hybrid shortest path and ray bending method for traveltime and raypath calculations[END_REF][START_REF] Zhao | Minimum travel time tree algorithm for seismic ray tracing: improvement in efficiency[END_REF][START_REF] Bai | 3D ray tracing using a modified shortest-path method[END_REF][START_REF] Bai | 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications[END_REF][START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF] (also called minimum travel time tree -MTTT) which are relatively underdeveloped algorithms based on Huygens' and Fermat's principles. Compared with the traditional ray tracing algorithms, the grid based scheme has several advantages for application to NDE modeling: it is able to easily locate ray paths in shadow zones; it is numerically stable and can find the first arrival in complex media.

The current study proposes, for modeling the propagation of an incident beam towards a receiving probe in a specimen of irregular surfaces including flaws, an approach based on seismic SPM algorithms adapted to ultrasonic NDE and specifically designed for the two points ray tracing problem. For a 2D propagation, the initial SPM Moser's algorithm [START_REF] Moser | Shortest path calculation of seismic rays[END_REF] consists in launching rays from a source, in considering in a 2D grid the adjacent node with minimal travel time as a secondary source, in repeating the previous process so as to propagate first arrival waves to the 2D grid nodes. More recently, Zhang et al [START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF] developed the interface source method (ISM): as in Moser's SPM, first arrival waves still propagate using a ray shooting scheme from nodes to nodes in a 2D grid but on the other hand only interface points are taken as secondary sources. The developed GIRT model uses a new SPM approach specifically devoted to two points ray tracing in a propagation medium constituted of homogeneous volumes (in terms of sound velocity). The GIRT takes advantages that in the latter case waves propagate in straight lines between interface points: only the specimen interfaces are meshed (in a simple 1D mesh for a 2D propagation), the source (on the emitter) and observation (on the receiver) points are added to the mesh and the wave ray path of minimum travel time passing through mesh points is found between the two fixed points. One of the major advantages of this GIRT method is the fast processing of very irregular surfaces as only the surfaces of the specimen are meshed and not all the specimen volume like in other SPM algorithms [START_REF] Nakanishi | A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure[END_REF][START_REF] Moser | Shortest path calculation of seismic rays[END_REF][START_REF] Cao | Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm[END_REF][START_REF] Van Avendonk | Hybrid shortest path and ray bending method for traveltime and raypath calculations[END_REF][START_REF] Zhao | Minimum travel time tree algorithm for seismic ray tracing: improvement in efficiency[END_REF][START_REF] Bai | 3D ray tracing using a modified shortest-path method[END_REF][START_REF] Bai | 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications[END_REF][START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF]. In addition, GIRT directly traces in an optimized manner (a two points ray tracing) the shortest path from a source reaching any required observation point located anywhere (outside the 2D grids needed by previous methods to propagate all shortest paths to all grid nodes). Furthermore, compared to the ISM method [START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF], the developed GIRT algorithm has notably added features to meet the specificities of ultrasonic NDT applications. Firstly, whereas ISM has been shown to deal with only the first arrival or one reflection path for P bulk waves, the GIRT algorithm is adapted to find not only the head wave ray path (first arrival) but also all complex later arrival waves including mode or nature conversion by adding constraints on the searched ray path, fixing one or more points through which the ray path must pass, and by defining independent types of wave propagation along each path connecting the constraint points. Secondly, the GIRT can take into account all types and natures of waves propagation occurring in NDE. Indeed the propagation is available not only for the SH head waves studied in [START_REF] Zhou | Ray path of head waves with irregular interfaces[END_REF] (also called TH) but also for P and SV waves (also called L and TV). Surface waves as Rayleigh waves and head waves are taken into account as well as diffractions from surfaces in any kind of waves. Finally, GIRT approach is extended to flaw scattering modeling since the flaws contained in the specimen are also meshed by the GIRT algorithm in order to model rays diffracted from flaws. Consequently the GIRT algorithm is generic for application to NDE as it can deal with any ultrasonic wave propagating near irregular surfaces or flaws of any geometry (CAD defined). Finally this tool is able to calculate valid geometric ray paths after interaction with interfaces or flaws. One interest of the study presented in this paper is also to evaluate the applicability of such a method initially based on seismic ray tracing in NDE applications for which the characteristic dimensions (propagation distance or scatterer size) to the wave length ratio can be smaller than in geophysics. The current study shows that GIRT leads to a better understanding of the head wave propagation and its interaction with the surfaces. In that aim, a data fusion of simulated results obtained by both GIRT and a FEM-based method is carried out.

The paper is organized as follows. The theoretical principles of the approach, the GIRT algorithm process and the GIRT capabilities are described in section 2; some ray tracings and wave fronts of the head wave and some later arrival waves, simulated with GIRT in 2D applications, are compared for validation with FEM results in section 3.

The Generic Interface Ray Tracing (GIRT) method

The development of a ray tracing method allows a better understanding of the head wave propagation on irregular surfaces by providing valid ray paths. The interest of such a study is described in parts 2.1 and 2.2. The developed ray tracing method (GIRT) is generic by the fact that it supports any kind of wave propagation and interaction with an irregular surface, and describes the physical phenomena involved in the head wave propagation. Its theoretical principles, functionalities and operating process are detailed in part 2.3.

Head waves in planar surfaces configurations

For a planar surface inspected in TOFD, a head wave is generated at the interface between the material and the emitter and propagates in the material along the surface, and is radiated as a bulk wave in the wedge, reaching the receiver. For planar surfaces, the head wave is seen as the result of a critical refraction phenomenon of the incident wave on the interface between the two media. Fig. 1 describes, in a TOFD inspection of a planar specimen containing an embedded crack, the propagation of different waves including the head wave. Fig. 2 illustrates the B-scan measured in the inspection described in Fig. 1. Three kinds of signals can be observed and are characteristics of a TOFD inspection on a planar specimen.

These signals correspond respectively to the head wave (1), the diffraction (2, 2') of the refracted bulk wave on the edges of the crack and the reflection (4) of the refracted bulk wave on the specimen bottom.

Head wave in irregular surfaces configurations

TOFD inspections of specimen with irregular entry surfaces can be required for specific NDE applications. For example, the process of repairing a material including a surface breaking crack implies to remove the matter around the defect then to replace the removed matter. The TOFD technique is applied on the specimen between the two reparation stages in order to ensure that the entire defect has been removed. The geometry of the specimen after matter removing is presented in Fig. 3.

A good understanding of the head wave propagation on such surfaces is essential for the complete modeling of TOFD inspections on realistic configurations. In particular, the head wave signal is one of the main references for the detection and the location of potential defects in the specimen. As explained in [START_REF] Zhou | Ray path of head waves with irregular interfaces[END_REF], irregularities of the specimen surface induce modifications on the head wave received signal compared to that observed on planar interfaces. Indeed such modifications can be noticed on experiments: the B-scan presented in Fig. 4 is extracted from a TOFD inspection measurement on the specimen described in Fig. 3 and highlights two kinds of differences which can be observed when scanning over the surface irregularities. Firstly, the head wave time of flight appears to be incoherent with the hypothesis of head wave propagation all along the entry surface: indeed the measured travel time is less than that of a surface wave propagating all along the complex surface. Secondly, the irregular surface causes an important decrease (more than 15dB) of the head wave amplitude compared to a planar interface.

To understand the physical process responsible for the observed head wave propagation on complex surfaces, an analysis of the nature of the field propagating near such surfaces has been carried out: the emitter generates a refracted field in the material which interacts with the surface. In particular, surface corners are responsible for diffracting the field in their geometrical shadow. Thus head wave propagation cannot be considered as a critical refraction anymore like in the case of planar surfaces, but the result of the diffractions of the generated refracted field on the surface irregularities. Fig 5a describes such behavior on an example of irregular surface: the head wave, after a bulk refraction of the incident wave, is diffracted two times on the surface corners, in order to create a diffracted wave, which is then refracted again towards the receiver. In order to visualize this behavior, some numerical simulations have been carried out with the hybrid technique CIVA/ATHENA [START_REF] Gengembre | A semi-analytic-FEM hybrid model for simulating UT configurations involving complicated interactions of waves with defects[END_REF] using both analytical models and finite elements simulations near the surface irregularities. The principle of this technique is the following. A volume in the specimen is defined, embedding the surface and its irregularities. The field reaching the volume entry is calculated by the analytical CIVA pencil method on the volume border. Then the field is propagated in the defined volume by the finite elements scheme ATHENA. The results of this simulation are presented in Fig. 5b and show that there exists a field diffracted from the surface corner in the geometrical shadow area of the irregular surface: the corresponding diffracted rays well generate the head wave (Fig 5a).

GIRT objectives and principles: construction of valid ray paths

Considering the behavior of the head wave near an irregular surface, the ray theory formalism is then applied for the two points ray tracing problem. The head wave propagation and more generally all the waves in the specimen are then the result of specific interactions between the wave field inside or outside the specimen and the specimen surface: the ray path representing this propagation is the combination of specific rays connected to each other by secondary source points. These secondary source points are located on the specimen surface and correspond to the locus of each interaction responsible of the wave propagation. The rays are a representation of the wave propagation between each secondary source points. 6

. The first configuration (Fig. 6a) is a three propagation media one for inspection of a specimen including a corner irregularity on its entry surface, an internal interface with slopes and a defect. The case of an isotropic elastic specimen with an entry cylindrical surface is shown in Fig. 6b. The upper medium is a fluid one in both cases. Each ray path in Fig. 6a and

i i M ∈ is consequently not exhaustive.
As to the configuration of Fig. 6a, the ray path Ray paths not shown on Fig. 6a and Fig. 6b result from other interactions in the specimen.

For example, mode conversion of the emitted L bulk wave occurs at the refraction on the entry surface, at specular reflections or at diffractions on the specimen surface irregularities or flaws (like at 8 M in Fig. 6a) and allows the propagation of several transversal bulk waves into the specimens. Mode conversions of bulk waves into Rayleigh-like waves also occurs on surface irregularities (at 4 M , 5 M , … in Fig. 6a or 3 M in Fig. 6b).

Finally the purpose of the GIRT is to find all physical paths described above in isotropic or anisotropic media using the Generalized Fermat's Principle. This principle is a basic one of an asymptotic ray theory, the Geometrical Theory of Diffraction (GTD) [START_REF] Keller | Geometrical theory of diffraction[END_REF] which extends Geometrical Acoustics by adding rays diffracted from canonical scatterers (edges, corners, curved smooth objects…) to the classical geometrical (refracted and reflected) rays. In the aim of ray tracing, the Generalized Fermat's Principle allows to determine if a ray path represents a physical wave propagation and is consequently valid. Indeed this principle indicates that the valid ray path for a wave resulting from particular interactions between the wave field and surfaces and is the one providing the stationary time of flight and respecting the predefined conditions imposed to the wave propagation (choice of the modes for instance). This principle is verified by finding the minimal time of flight and the associated ray path. Consequently the principle of the GIRT is to find the valid ray path respecting the Generalized Fermat's principle and achieving the wave propagation conditions using a time of flight minimization technique. Its algorithm process is described in the following section.

GIRT algorithm process

The steps of GIRT algorithm process are described in this section for the general case of anisotropic media and schematically represented in the functional diagram of Fig. 7. This algorithm process is the subject of a patent filing [38].

This process is first described when used for determination of the head wave path, the first arrival (like the paths 1 on Fig. 6a and Fig. 6b) and is called in that case GIRT-A. The input data (step 1) of the GIRT are the CAO description of the specimen in which the ray path has to be calculated, the source and observation points which are the ray extremity points, the slowness surfaces for the wave propagation mode desired inside and outside the specimen, the description of defects inside the specimen and the discretization fineness of the specimen or flaw surfaces which defines the calculation accuracy.

Indeed the next steps of the algorithm process is the discretization of both each specimen surface (step 2) and of each existing defect (step 3) in order to deal with any irregular surface:

for the configuration of Fig. 6a, the two irregular surfaces and the planar defect are meshed.

The obtained discretization points collection { } (1,..,N) i i M ∈ constitutes the secondary sources collection for waves diffraction as stated by the Huygens principle in a propagation medium constituted of homogeneous volumes. Consequently these secondary sources are possible i M points composing the ray path as described in section 2.3. The GIRT at this step 2 uses a similar approach to the interface source method (ISM) developed by Zhang and al [START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF] since interface points are seen as secondary sources.. On the other hand, the discretization carried out by GIRT differs from that performed by ISM [START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF] and other SPM algorithms [START_REF] Moser | Shortest path calculation of seismic rays[END_REF]. Indeed as only the surfaces and the defects are discretized, the ray path searching is done only through secondary sources located along the specimen and defects surfaces and not trough a grid covering all the specimen volume. The obtained secondary sources collection has a limited size N compared to that generated from a volume discretization and provide a more Method (SPM) [START_REF] Moser | Shortest path calculation of seismic rays[END_REF] in the seismic domain and is applied in the GIRT (step 7). The secondary source points collection are the nodes of an oriented graph. In this graph, two different vertices are connected by an elementary ray path included in the collection. At each path is associated a cost which is the elementary time of flight between of the path. Once the graph is obtained, the optimal route (leading to the shortest time of flight between the source and observation points) is obtained by the classical Dijkstra algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF], which has the major advantage to be a fast route search algorithm since all possible routes are not calculated: only the most likely optimal routes are taken into account. At that step, the GIRT algorithm differs strongly from the ISM developed by Zhang and al [START_REF] Zhang | A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources[END_REF]. ISM is conceived to perform ray shooting in a grid meshing all the volume by launching rays from a source in all directions and by propagating first arrival waves from nodes to nodes in the grid: consequently a ray path reaching a predefined observation point cannot be found directly notably if this point is not a grid node. In GIRT, the application of the Dijkstra algorithm on the secondary sources collection is well fitted and specific of the two point ray tracing problem. Indeed, the source and observation points have been integrated in the secondary sources collection so that the fast route search gives directly the valid ray path between these two points.

Later arrival waves (for example all paths represented in Fig. 6a and Fig. 6b except the paths 1) can also be modeled by a GIRT-B version by adding constraints on the searched ray path, fixing one or more points through which the ray path must pass. In the case of a ray path calculation with obligatory points of passage, the GIRT-B process differs from that previously described for the head wave ray path (GIRT-A process) only by some slight changes. First, the GIRT-B inputs include in addition the list of sorted constraint points and the definition of independent types of wave propagation along each path connecting these points. Then the GIRT-B operates as follows: a couple of points is formed by the source point and the first constraint point, then several other couples by two adjacent constraint points and the last couple by the last constraint point and the observation point. An optimized ray is then searched between the two points of each couple: in that aim, the previous points are seen as the source and observation points for GIRT-A process and the GIRT-A process is consequently carried out independently for each couple, the latter having its own type of propagation, allowing all type of wave propagation in the specimen. At last, all the calculated couples ray paths are concatenated to obtain the complete ray path between the source and observation points.

GIRT capabilities

GIRT is thus a generic algorithm calculating geometrical paths and the corresponding ray time of flight propagating in isotropic or anisotropic media. It has been designed to determine the valid ray path between two points assuming one or several given interactions between the current ray and the irregular surface. The two points can be located anywhere inside or outside the specimen.

Two kinds of interactions, accounted by GTD [START_REF] Keller | Geometrical theory of diffraction[END_REF], are supported by GIRT: diffraction on smooth objects, leading to surface creeping rays (see the ray 5 6 M M on Fig. 6b), and diffraction on wedges or edges (see Fig. 6a), leading to bulk or surface diffracted rays.

Furthermore, the propagation of numerous wave modes is taken into account with potential mode conversion at the surface: the longitudinal and transversal modes (possibly occurring for different kinds of waves: bulk waves, lateral waves, reradiated head waves or surface waves…), and the surface Rayleigh-like waves.

All the ray paths presented on Fig. 6a and Fig. 6b are modeled by setting some wave propagation conditions supported by the GIRT. As an example the head wave ray paths 1 of

GIRT simulated results and theoretical validations

Calculations using GIRT of ray paths and wave fronts are shown in this section in the aim of modeling the propagation of head waves near irregular surfaces in some TOFD configurations. Such simulations are compared and superimposed to FEM ones.

Validation of the head wave time of flight

A theoretical validation of GIRT is performed by comparing times of flight calculated with both GIRT and finite-elements (FEM) simulations obtained with CIVA/Athena, which is considered here as a reference model. The comparison deals with TOFD inspection simulations for several inspection directions (θ = 45°, 60° and 75° where θ is the angle defined in Fig. 8) with L waves emitted at 5MHz. The specimen is a planar stainless steel block which is inspected using two contact planar rectangular transducers (6mm*5mm). The specimen entry surface has two slopes. The head wave path determined by GIRT is illustrated in blue in Fig. 8 and includes a diffraction on the surface corner.

The times of flight obtained by GIRT have been compared to the results of the FEM simulation software for different slope angles and are shown in Fig. 9. Each graph of the Fig. 9 is associated to a particular inspection direction (from L45° to L75°), the solid lines are the results obtained with the GIRT and the dashed lines are the results obtained by FEM simulations.

The graphs show a very good concordance between the results of the GIRT and FEM simulation, for each direction of inspection and any slope angle. The GIRT is then able to calculate the head wave time of flight with a very good precision (error less than 0.1%). The simulated times of flight of the head wave signal received by the probe are almost independent of the inspection directions because the head wave signal is generated by only some refracted rays included in the beam radiated in the specimen by the emitter surface, those rays being close to the mean head wave path shown in blue in Fig. 8 including a diffraction on the corner.

Analysis of the propagation of waves in the surface shadow

Additional finite-element simulations on CIVA/Athena have been performed on irregular surfaces including curved parts, and the propagation wave field in the material at different times is also visualized on FEM snapshots. In this section, GIRT has been employed in the aim of identifying the nature of waves propagation corresponding to some observed wavefronts in the FEM snapshots.

Indeed, when studying the interaction on a beam with an irregular surface, numerous wave fronts are generated. For instance Fig. 10 shows FEM simulation of the waves propagation near a curved valley constituted of two curved edges separated by a planar part. Numerous wave fronts are observed in the FEM snapshot of Fig. 10, which highlights the great number of interactions occurring between the refracted field and the specimen. The nature of some of these wave fronts (numerated in Fig. 10) can be easily guessed. The surrounded front L1 corresponds to the direct L-wave refracted in the specimen by the emitter. The front L1 is reflected by the specimen bottom (dashed black horizontal line in Fig. 10) creating the Lwave front L2 and T-wave front T1. T2 refers to T-waves refracted in the specimen. Other wave fronts are of more complex nature. It is consequently of interest to identify the nature of rays which have generated the observed wave fronts in the FEM simulations.

So as to identify the nature of one particular FEM wave front, the following methodology is chosen:

-The studied wave front is identified on the FEM snapshot as the local maximum of the displacement field.

-An appropriate hypothesis on the wave propagation and interactions with the surface is formulated by the GIRT user to explain the generation of the studied wave front.

Considering the previously chosen path hypothesis, GIRT builds rays until their time of flight reach the FEM snapshot time.

-The resulting wave front simulated by GIRT is then drawn for the FEM snapshot time by tracing the curve orthogonal to every ray found at this time of flight [START_REF] Bouche | Asymptotic methods in electromagnetics[END_REF]. The GIRT simulated wave front is superimposed to the FEM snapshot and its location compared to that of the FEM studied wave front to validate the chosen propagation hypothesis.

In the present study, identification of wave fronts nature by comparing FEM and GIRT simulations is done only for wave fronts relevant for studying head waves propagation. But a similar analysis could be performed for all the other wave fronts observed on FEM snapshots.

The method is thus applied for the following wave fronts:

-The front of the refracted L-wave (Fig. 10), whose diffraction on the second curved part of the surface specimen is responsible of head wave propagation towards the receiver (as said in section 2.2).

-The L head wave (Fig. 11 and13), received at first on the receiver.

-The T wave received just after the L head wave on the receiver.

For each figure, the ray paths (white lines) and the front (dashed black line) simulated by GIRT are superimposed to the studied FEM snapshot (in color code).

In (at the interface specimen/wedge of the receiver) and give the results of respective Fig. 11 and in Fig. 12. In Fig. 11, a second refraction of longitudinal waves without mode conversion is assumed and GIRT-A is used to calculate the shortest rays between the source an observation points located on the receiver crystal surface. In Fig. 12, a second refraction with conversion of the L wave into a T wave is supposed. To apply the last hypothesis, a constraint point is defined on the plane part of the irregularity and the GIRT-B algorithm is used so as to calculate the rays using the following assumptions: when the path takes place inside the specimen, it is always travelled at the celerity of the L wave; when the path takes place both outside the specimen and from the source point to the constraint point (i.e. in the emitter wedge), it is travelled at the celerity of the L wave; when the path takes place both outside the specimen and from the constraint point to the receiver point (i.e. in the receiver wedge), it is travelled at the celerity of the T wave. By this way ray paths calculated by the GIRT-B satisfy the mode conversion L-T condition on the surface near the receiver. As in Fig. 10 and for the two cases in Fig. 11 and Fig. 12, wave fronts calculated by GIRT perfectly match with two different wave fronts simulated by FEM. Furthermore the wave of type T (Fig. 12) is received by the receiver later than that of type L (Fig. To emphasize this conclusion about the nature of the head wave propagation, the GIRT simulated rays of a L-surface wave propagation all along the irregular specimen surface without any specimen bulk propagation hypothesis are shown in Fig. 13. To satisfy the previous condition, two constraints points are specified to the GIRT-B algorithm and are positioned on the specimen surface respectively just before and just after the surface irregularity, so as to impose the propagation of rays along the irregular surface. In that case, the surface ray is generated by critical incidence on a planar interface as described in Fig. 1.

The corresponding extracted wave front does not match the head wave front obtained from FEM simulations and confirms that the bulk mechanism stated previously for the head wave propagation is likely to predominate for very corrugated surfaces.

The results in the previous comparisons between GIRT and FEM show that the head wave propagation in the bulk implied by complex interactions between the wave and the surface must be taken into account to correctly model the head wave propagation on irregular surfaces. These comparisons between the developed ray tracing model (GIRT) and FEM simulations validate the choice of the GTD approach used in GIRT to model the complex interactions between an irregular surface and the waves radiated by the probe. The GIRT based on GTD principles is able to compute valid ray paths providing wave fronts in agreement with FEM simulated results in terms of both location and shape.

Conclusion

The modeling of TOFD configurations used in ultrasonic NDT requires the study of the head wave propagation. While analytical theories for head waves on planar interfaces are well established, the complete simulation of head waves on irregular surfaces has to be developed.

In that aim, one suitable approach is to use the Geometrical Theory of Diffraction.

A generic algorithm of ray tracing between interface points (GIRT) has thus been developed using one GTD principle: the Generalized Fermat's Principle. GIRT is a specific and generic adaptation of seismic ray tracings to NDE modeling of waves propagation and to the two points ray tracing problem solving: a fast processing of the waves propagation near very irregular surfaces between two fixed points is obtained by meshing only those surfaces and the defects rather than all the volume of the propagation media like in existing seismic ray tracings. The approach is adequate for a propagation medium constituted of homogenous or weakly inhomogeneous volumes: in the latter case, the wave propagation can be modeled by tracing straight lines in a mean homogenous volume and adding travel time fluctuations [START_REF] Lu | Stochastic simulation of the high-frequency wave propagation in a random medium[END_REF].

The advantages of the GIRT are the specific treatment of diffracted and creeping rays in shadow areas, the account of both any irregular surfaces of the specimen and embedded defects (both CAD defined), the wave mode conversions at interfaces or flaws, and the calculation of all waves (bulk, surface and head waves) propagating towards a predefined observation point during a TOFD inspection for NDE including both first arrival (head wave) and latter arrival waves.

By modeling the complex wave interactions with the surfaces thanks to this ray approach, GIRT is able to correctly predict the head wave time of flight and to interpret the complex waves propagation at the vicinity of an irregular interface. For very corrugated surfaces, a bulk propagation mechanism has to be taken into account for the head wave modeling.

Valid geometrical rays paths are henceforth determined for head waves using GIRT pure kinematic analysis. Use of GTD for simulation of the head wave amplitude will be investigated thereafter by modeling surface creeping rays and corner diffractions. GIRT could then lead in the future to a complete modeling of the head wave propagation on irregular interfaces.

Figure 4 (b&w)

Figure 5 (b&w)

Figure 6 (b&w)

Figure 7 (b&w)

Figure 8 (b&w) • This article provides a modeling study of head waves near irregular surfaces in NDT.

• Head wave propagation near such complex surfaces implies bulk mechanisms.

• A generic algorithm of ray tracing between interface points (GIRT) is developed.

• GIRT, based on Generalized Fermat's Principle, models effects of complex surfaces.

• The head wave fronts computed by GIRT are in good agreement with FEM simulations.

2 (∈

 2 The GIRT ray tracing formalism is thus based on the Huygens principle. Considering known source point s M and observation point o M , there exists a secondary source points collection { } (1,..,N) i i M ∈ located on the specimen/flaws surfaces which represents the locus of all potential interactions between the wave field and the surfaces: the size N of the secondary source points increases with the complexity of the surface irregularities. The ray path associated to a wave propagation between the points s M and o M is then a combination of elementary particular rays ( ) 2 ( , ) (s,o,1,.., ) j l j l N M M ∈ whose extremities are secondary source points i M among the collection { } (1,..,N) i i M ∈ of an elementary ray emitted from s M and of an elementary received on o M . In the following, it is illustrated on examples (Fig. 6) that this principle of secondary sources can be applied to represent all physical waves in a propagation medium containing interfaces and flaws. In a general manner, numerous waves including the head wave are emitted from a source point s M and received on an observation point o M : these waves are characterized by their interactions with the specimen surfaces involved in their propagation and then by the secondary source points constructing the associated ray path. Moreover, theses interactions may be responsible of a wave mode or wave nature conversion: consequently, each elementary ray ( ) can represent the propagation of a longitudinal (L) bulk wave, transversal bulk waves (T in an isotropic medium, qT1 and qT2 in anisotropic one), lateral surface wave (L or T) or Rayleigh-like (Rayleigh, Generalized Rayleigh, Stoneley) wave depending on the type of interaction located on the secondary source point j M which is the departure point of ray ( ) 2 ( , ) (s,o,1,.., ) j l j l N M M ∈ . Two inspection configurations involving irregular surfaces are studied as examples in Fig.

Fig. 6b represents

  Fig.6brepresents one particular wave propagating. Only some ray paths among all existing

  M M (path 1) is the ray with the shortest time of flight and thus corresponds to the head wave propagation: its path is composed of a refraction at 2 M of the longitudinal bulk wave emitted at s M , a corner diffraction at 4 M without mode or nature conversion and a refraction towards the observation point o M . However the diffraction of the L bulk wave at the corner 5 M implies also a mode conversion into a T wave and gives after refraction at 11 M M M (path 2 in Fig. 6a) with a higher time of flight than the head wave ray path. The ray path M M M M (path 3) is associated to a L bulk wave critically refracted at 3 M into a L lateral wave which propagates along the surface and is critically reradiated at 10 M into a bulk L wave. As another example the secondary source point 7 M also converts the lateral wave into a Rayleigh wave which radiates at 9 M a bulk L wave and gives the ray path M M M (path 5) corresponds to a L bulk wave emitted at s M and refracted at 1 M which is reflected on the internal interface of the specimen on 6 M and diffracted on the edge defect on 8 M before reaching o M .

Fig. 6b shows

  Fig.6bshows the ray paths of different waves which propagate close to a surface with a

2 (

 2 stable and faster research of ray paths.To determine the valid secondary sources involved in the searched ray path in the view of the two points ray tracing problem, the elementary paths between two secondary sources are calculated (step 4): an elementary path represents the wave propagation along a direct path involving two and only two i M secondary sources without any other diffraction effect. All the valid elementary paths between two i M points respecting the previous condition are calculated for all the i M points and the associated elementary times of flight are deducted from the slowness surfaces given and stored by the algorithm process (step 5). As an example, the ray 2 5 M M in Fig.6bis an elementary ray path of the path 1 and path 2. A secondary sources collection, an elementary paths collection and an associated elementary times of flight collection have been obtained during the previous steps. The s M source and o M observation points are added in step 6 to the existing secondary sources collection and all the new valid elementary paths involving one of these two points are calculated and taken into account in the collections. The interest of having separated this step 6 from step 5 is to allow the ray path calculation for different couples of source and observation points but with the same specimen without needing to recalculate at each ray calculation all the other elementary paths in the specimen (already independently determined at step 5). Finally, knowing the secondary sources collection and the elementary paths collection, the search for the valid wave ray path in the view of the two point ray tracing reduces to an optimization problem. It indeed corresponds to the determination of which secondary sources i M are involved in the ray path and then which elementary rays ( ) associated to the minimal time of flight. This approach is called the Shortest Path

Fig. 6a and

 6a Fig.6aand Fig.6bare calculated using by finding the shortest ray path for a longitudinal

Fig. 10 ,

 10 ray paths simulated by GIRT are obtained considering the hypothesis injected in GIRT of a diffraction of the refracted L-wave field on a curved object without mode conversion and without constraint point. Ray paths are thus built with the GIRT-A algorithm between one source point (emitter center) and a set of observation points located in the specimen bulk after the second curved part at right, inside the shadow areas formed by the surface irregularity. The GIRT objective is here to determine between the source and each observation point the shortest path which corresponds to the head wave path. The calculated ray path indicates that the emitted wave is refracted in the specimen bulk, before being diffracted by the first curved object at left. The GIRT modeling of the diffracted field by near this curved part reveals creeping rays. Indeed, the refracted longitudinal rays reaching tangentially the first curved object are converted into creeping rays. These creeping rays are then diffracted into both bulk rays along their curved propagation and a surface ray which propagates along the planar surface. The surface ray is converted again at the junction with the second curved part into a creeping ray which is diffracted into bulk rays generated tangentially to the second curved object at right inside the shadow area caused by the surface irregularity. As seen in Fig.10, the calculated wave front is the result of the wave diffractions on the two curved objects. The previously described head wave propagation is validated by comparing the GIRT and FEM simulated wave fronts and observing a good agreement between them.The following two cases, presented inFig 11 and Fig 12 are obtained with the same configuration as in Fig. 10, but the propagation of waves is extended up to the receiver, allowing observing the wave fronts received in experimental studies. As in Fig. 10, the refracted longitudinal wave is diffracted by the curved parts at left and at right and refracted again at the last surface between the specimen and the wedge of the receiving probe. Two different hypotheses detailed hereafter are inputted in GIRT concerning the second refraction

11 )

 11 as shown by comparing the location of their corresponding wave fronts for a fixed time of flight. Consequently, in the previous configuration, the head wave observed at the receiver which is by definition the first received wave logically mainly corresponds to a complete path without mode conversion involving only faster L waves. All the other FEM simulated wave fronts obtained in the snapshot of Fig. 11 lead to later arrival times if they reach the receiver.

Figure 1 :Figure 2 :

 12 Figure 1: A NDT inspection using TOFD technique

Figure 3 :Figure 4 :

 34 Figure 3: An irregular geometry studied by TOFD technique

Figure 5 :

 5 Figure 5: Propagation of a HW on an irregular surface (valley consisting in several planar

Figure 6 :

 6 Figure 6: Various paths (color lines) resulting from several interactions with two kinds of

Figure 7 :Figure 8 :

 78 Figure 7: Functional diagram of the GIRT algorithm process

Figure 9 :

 9 Figure 9: Comparison of times of flight calculated both with GIRT and with FEM simulations

Figure 10 :

 10 Figure 10: Diffraction of L-waves on curved objects. Head wave ray paths calculated by

Figure 11 :

 11 Figure 11: Head wave (L wave) received by the transducer. At right: zoom on the receiver.

Figure 12 :

 12 Figure 12: Hypothesis of a T wave received by the transducer. At right: zoom on the receiver.

Figure 13 :

 13 Figure 13: Hypothesis of a L wave received by the transducer without propagation in the

Figure 1 Figure 1 Figure 2 Figure 2 Figure 3 Figure 3

 112233 Figure 1

Figure 9

 9 Figure 9 (b&w)


Romain Coulon

Vladimir Kondrasovs

Karim Boudergui

Stéphane Normand

Moving Sources Detection Algorithm for Radiation Portal Monitors Used in a Linear Network

Keywords: Algorithm, detection, filter, instrument, measurement, monitor, moving, nuclear, radiation, radioactive, safeguard, security, signal processing, source, time series analysis

To monitor radioactivity passing through a vehicle such as a pedestrian, a car, a train or a truck, Radiation Portal Monitors (RMP) are commonly employed. These detection systems consist of a large volume detector set close to the potential source path. An alarm is then triggered when the signal rises over a threshold initially estimated in view of the natural background signal. The approach developed in this work makes use of several detectors in a network along the source path. The correlation detection approach is elaborated to take into account the temporal periodicity of the signals taken by all distributed sensors as a whole. This new detection method is then not based only on counting statistics but also on the temporal series analysis. Therefore, a specific algorithm has been developed in our laboratory for this security application and shows a significant improvement, especially in terms of detection probability increase and false alarm reduction. This paper presents the theoretical approach and promising results obtained by simulation.

I. INTRODUCTION

M ONITORING a moving source is currently required for numerous nuclear security applications. Nuclear materials could be contained into a pedestrian, a car, a truck, or a t r a i na n dc o u l da l s ob ed i s s e m i n a t e di n t oaw a t e rp i p eo ra n air vent. The detection is challenging due to the source motion when passing by the probe.

A. State of the Art

To optimize the detection eff ciency, a large volume sensor is required. Such systems are proposed by a large number of industrial suppliers [START_REF]Radiological control of vehicle load[END_REF], [START_REF]RADOS DnD-Portal[END_REF]. Advanced RPMs, additionally, are currently developed in research laboratories, exploiting the information contained into the deposited energy distribution given by the large plastic scintillator [START_REF] Kangas | The use of artif cial neural networks in PVT-based radiation portal monitors[END_REF]- [START_REF] Runkle | Analysis of spectroscopic radiation portal monitor data using principal components analysis[END_REF]. Time series analysis of the counting signal is also studied to improve the detection of a single RMP [START_REF] Robinson | Time series evaluation of radiation portal monitor data for point source detection[END_REF], [START_REF] Vilim | Sensitivity improvement in low-prof le distributed detector systems for tracking sources in transit[END_REF]. The work developed in this paper deals with time series analysis applied to RPMs in linear network.

The authors are with the CEA, LIST, Laboratoire Capteurs et Architectures Électroniques and Laboratoire, F-91191 Gif-sur-Yvette, France (e-mail: romain.coulon@cea.fr; vladimir.kondrasovs@cea.fr; karim.boudergui@cea.fr; stephane.normand@cea.fr).

B. Conventional Approach

Conventional RPMs are implemented with a single channel and gross count rate continuous measurements. Detection is commonly done from a given channel by a threshold. The threshold is calculated according to the background f uctuation band on the backed level estimated with a Poisson counting assumption and an acceptable false alarm probability FA P . Counts from background are considered as a Poisson process. Its expectation and its variance are a function of background eff ciency , background intensity around the detector and elementary counting time as described in Eq. ( 1).

(

) 1
An decision threshold DT is determined by the statistical standard deviation and a coverage factor (set in regards of the maximum acceptable FA P ).

(

Counts from signal are similarly considered as a Poisson process. Its expectation and its variance are the functions of the source eff ciency , the source activity and elementary counting measurement time as described in Eq. ( 3).

(

) 3
The detection test consists in trigging the alarm if the condition becomes true. The ratio is therefore proportional to the source eff ciency and activity, proportional to the square root of the elementary counting time and inversely proportional to the square root of the background eff ciency and intensity as shown in the Eq. ( 4).

(4)

The background eff ciency increases proportionally to detector volume. The source eff ciency converges to an asymptotic value as the detector volume increases. An optimum in terms of signal to noise ratio is calculated (by MC particles transport code for example) as a function of the experimental conditions (nature of the source, background, vehicle…). Consequently, optimizing the design of the sensor comes f rst in the issue to be addressed.

The counting time has to be short compared to the speed of the source. In practice, it can only by optimized by making the considered vehicle stop in front of the system. Such a method could disturb the traff c f ux; it is therefore undesirable.

C. Concept of the Correlation Approach

Another way of improvement is proposed in this work: it relies upon several portal units set along the source path. As a result, the detection probability is increased and the network structure could provide us with additional type of information about the signal thanks to the temporal correlation between sensors. Instead of trigging on the independent signal amplitude, a periodicity into the overall signal is searched for. An improvement in detection ability is suggested. Moreover, this new approach could be more reliable as far as the self-shielding issues are concerned.

Past studies have been carried out for the source tracking using matrix sensor networks. MCNC calculations and Bayesian approach allowing sources tracking without any knowledge of the source trajectory have been developed [START_REF] Brennan | Radioactive source detect i o nb ys e n s o rn e t w r k s[END_REF]. An algorithm based on triangulation has been studied for source localization purposes [START_REF] Chin | Identif cation of low-level point radioactive source using a sensor network[END_REF]. These methods are optimized for 2-D problems and are therefore unsuited for our situation where the source trajectory is well known (1-D). Faster algorithms have then been developed but these require an accurate knowledge of the source velocity which is not usually the case for security problematic [START_REF] Nemzek | Distributed sensor networks for detection of mobile radioactive sources[END_REF], [START_REF] Stephens | Detection of moving radioactive sources using sensor networks[END_REF]. Fusion of correlated decisions is another path to address this issue. Local decision from multiples sensors are sending to a fusion center and correlation test is done (CUMSUM test or test based on copula theory) [START_REF] Quian | Quickest detection of nuclear radiation using a sensor network[END_REF], [START_REF] Sundaresan | Distributed detection of a nuclear radioactive source using fusion of correlated decisions[END_REF].

The next parts of this paper present a detection method based on temporal series analysis associated with linear sensor network without taking in account the local decision information. It allows a continuous and fast detection without accurate knowledge of the source velocity. The method is described and a comparative study with conventional detection method is presented.

II. METHODS

The radioactive source is passed in front of the N successive channels. Counts are recorded to a memory corresponding to every channel (see Fig. 1).

The detection f lter and measurement are processed with the following steps.

A. Memory Management

At the end of the counting time , a "First In First Out" (FIFO) list is recorded at every channel. After, the last counting is recorded into the f rst values of the memories. Counting values are noted where is the channel number and the time step at the memory level (see Fig. 2).

B. Correlation Vector Calculation

A temporal correlation vector is calculated as a function of delay times between sensors. The explicit formula is presented in Eq. ( 5).

(

With and .

C. Determination of the Decision Threshold

A decision threshold is calculated using the mean of the static product (Eq. ( 6)) and its associated empirical standard deviation (Eq. ( 7)).

The threshold estimation could be impacted by the source signal itself. To avoid this phenomenon the median could be used instead of the mean for the static product estimation. Another empirical standard deviation is recalculated after rejection of extreme values has been performed.

D. Detection Test

A detection test consists in determining if a correlation between the channels is considered relevant. The maximum of correlation is searched into the vector andatestisperformedby comparison of this value to the static product and its standard deviation. The parameter is the coverage factor set in regards of the acceptable FA P . This parameter can be considered as the threshold set proportionally the f uctuation level of the correlation product. [START_REF] Brennan | Radioactive source detect i o nb ys e n s o rn e t w r k s[END_REF] The system is able to detect all source motions as long as the speed can be considered constant and comprise into the delay time scan range. The system is also able to assure a reliable detection behavior during the evolution of the background signal. This approach is therefore immune to self-shielding effect.

E. Speed Estimation

If detection has occurred, the time period corresponding to the maximum of vector could give, using the distance between sensors , an approximation of the source speed .The accuracy in this speed estimation is linked to the time sampling of the detection system.

(

) 9
(10)

F. Signal Simulation

The signal is simulated by the superimposition of a constant background intensity and a source intensity (equal to the count rate in front of a detector). A geometrical model allows for the simulation of the global intensity of a channel as a function of time (see Fig. 3).

The intensity is derived as a function of the distance between probes , the distance to the source track line and the source velocity (see Eq. ( 11)-( 12)).

(11) [START_REF] Quian | Quickest detection of nuclear radiation using a sensor network[END_REF] The inhomogeneous Poisson process is generated for each channel as illustrated in Fig. 4. Pulse trains are generated using the pseudo-random function rand of Matlab into the exponential probability function modeling the pulse-time separation as a function of the signal intensity . The counts are integrated during an elementary counting time to obtain three independent count signals as presented in Fig. 5.

The signal is then characterized by the background intensity , the source intensity , the source velocity , the distance between the source and the detector line path ,andthecorresponding number of samples separating detector to each other's ,a s .

III. RESULTS

A. Detection Filter Testing

The correlation production is then calculated for several delay times as described in Eq. ( 5) and shown in the Fig. 6. The detection threshold (dash line) is calculated in regards to the signal correlation product f uctuations (see Eq. ( 6)-( 8)). In comparison to Fig. 5 (same displayed signals) obvious improvements in the signal to noise ratio are noticed when using the correlation product estimator.

To evaluate the performance of this scheme, the following detection algorithm conf gurations are simulated:

•T h e f rst one is the conventional method described in Eq. ( 1)-( 3) for a sensor network from 1 to 5 units. Detection occurs if at least one of the signals falls above the threshold. • The second one is the time correlation algorithm as described in Eq. ( 5)-( 8) with a sensor network of 1 to 5 units.

B. First Results

To evaluate FA P during the potential passage of the source, constant signals are injected ten thousand times during an interval of (corresponding to the velocity of the potential source used for DP calculation). This sequence is then generated times with an elementary sampling time of 0.1 second. The FA P is estimated by the number of detection occurring divided to the number of generated sequences (Eq. ( 14)). The associated statistical error is shown in Eq. ( 15). ( 14) [START_REF]EJ-200 plastic scintillator[END_REF] The FA P evolution as a function of coverage factor is shown in Fig. 7 for the two detection f lters.

The FA P distribution is completely different between the individual and the correlation detection methods. Signal from the source ( cps, cps, m.s m) is then injected ten thousand times to estimate the detection probability DP (see Eq. ( 16), ( 17)). ( 16) [START_REF] Kondrasovs | An adaptive smoother for counting measurements[END_REF] Fig. 8 shows the DP as a function of the coverage factor for both detection f lters.

The DP distribution for the correlation f lter is smoother than the statistical one. The detection probability from source SDP is def ned as the difference between the DP and the FA P ,which is a satisfactory indicator performance (see Eq. ( 18), (19)).

(

) (19) 18
The SDP evolution as a function of the coverage factor is presented in Fig. 9.

It can be seen that the optimal coverage factor is different between the two detection methods and in regards to the number of used channels. In this example, the correlation method clearly shows a signif cant detection probability improvement in comparison with the individual statistical threshold method.

C. Impact of Parameters

When an increase of the radioactive source speed occurs, the SDP decreases in the same manner for both detection f lter methods (see Fig. 10).

The distance between detectors does not impact the individual detection but could impact the correlation detection. One must therefore ensure that this distance is long enough within the speed range at hand.

D. Realization

A prototype system of the RPMs network is under development in our laboratory in the context of the SECUR-ED European project for urban mass transportation [START_REF] Boudergui | New monitoring system to detect a radioactive material in motion[END_REF].

Each sensor unit is composed of two plastic scintillation detectors coupled to fast linear focused 12 stages PMT. Scintillators are EJ-200 produced by Eljen Technology with the dimensions mm. Data acquisition electronics (DAQ) have been developed by CEA and an embedded video camera to record movies of the suspect object has also been implemented [START_REF]EJ-200 plastic scintillator[END_REF], [START_REF] Carasco | D a t a acquisition and analysis of the UNCOSS underwater explosive neutron sensor[END_REF].

The DAQ is made of a compact set including a DC module providing power supply to the other modules (5 W), a Field Programmable Gate Array (FPGA) board and a High Voltage board for the 2 PMTs.

The FPGA card aims at digitalizing the signal directly from the detector anode signal. Four shifted ADC with Msample.s are used to reach the 800 MHz sampling frequency.

T h es o f t w a r eu s e di nt h i ss y s t e m is based on a client/server model. The server software aims at managing and collecting data from sensor units, saving them into a f le and sending them to the connected client.

The client can be a portable computer (PC) based near the control unit or deported to a control room. A standard network is used to establish the connection and to communicate.

The server can collect all data every 100 ms from every sensor unit, perform the signal processing and send automatically an alarm in case of the radioactive source detection.

IV. DISCUSSION

The correlation of network portal monitor signals forms a promising approach for moving source detection applications. The distance between probes has to be carefully set by the user into the network. Adding a smoothering stage is a path for improvement. Future work will therefore approach by coupling time analysis detection with an upstream adaptive smoother [START_REF] Kondrasovs | An adaptive smoother for counting measurements[END_REF].

Because they do not scan the same type of information, the conventional and the correlation methods are not antagonistic but complement each other. We shall therefore conduct a study about the combination of both approaches.

Another way of improvement lies within the analysis of the evolution in the temporal prof le of the measurement period. The exploitation of this time delay evolution prof le is the next step in detection enhancement allowed by correlation approach.

V. C ONCLUSION

The detection of moving radioactive sources constitutes a challenge for nuclear security applications. Processing of the signal given by a RPM network has been studied. An algorithm allowing detection based on the time correlation between channels has been developed. A signif cant detection probability improvement has been shown by simulation study of RPM signals.

This preliminary result proves that time correlation approach is an innovative solution for moving source detection. Further improvements are envisaged on the basis of this algorithm, notably by coupling with an upstream adaptive smoother, coupling with the conventional individual f lter, and time delay evolution prof le approach. This algorithm will be implemented into a prototype and experimentally tested. A demonstration will be conducted during the SECUR-ED project demonstration, in order to verify compliance with international standards ANSI N42.35 [START_REF]American national standard for evaluation and performance of radiation detection portal monitors for use in homeland security[END_REF].

Fig. 1 .

 1 Fig. 1. Principle of the signal recording.

Fig. 2 .

 2 Fig. 2. Formulation of the signal matrix.

Fig. 3 .

 3 Fig. 3. Geometrical model.

Fig. 4 .

 4 Fig. 4. Example of pulses arrivals from 3 sensors while a source is passing.

Fig. 5 .

 5 Fig. 5. Simulated counts signal from 3 sensors while a source is passing.

Fig. 6 .

 6 Fig. 6. Correlation production as a function of the delay time during the pass of a radioactive source.

Fig. 7 .

 7 Fig. 7. FA P as a function of coverage factor for the individual detection f lter (double line) and the correlation detection f lter (single line).

Fig. 8 .

 8 Fig. 8. DP as a function of coverage factor for the individual detection f lter (double line) and the correlation detection f lter (single line).

Fig. 9 .

 9 Fig. 9. SDP as a function of coverage factor for a signal of 12 cps and a baseline background of 20 cps.

Fig. 10 .

 10 Fig. 10. SDP as a function of the source speed for a signal of 12 cps, a baseline background of 20 cps, distance between detector of 3 m, a coverage factor of 3 for the individual detection, and a coverage factor of 7 for the correlation detection.

Fig. 11 .

 11 Fig. 11. RPMs network prototype showing three long detection units and a small control unit.

ACKNOWLEDGMENT

The authors would like to thank J, Dumazert for his help in writing this paper.


Cindy Le Loirec
email: cindy.le-loirec@cea.fr

Dominique Chambellan

David Tisseur

IMAGE-GUIDED TREATMENT USING AN X-RAY THERAPY UNIT AND GOLD NANOPARTICLES: TEST OF CONCEPT

Gold nanoparticles (GNPs) have the potential to enhance the radiation dose locally in conjunction with kV X-rays used for radiation therapy. As for other radiotherapy modalities, the absorbed dose needs to be controlled. To do that, it is an advantage to know the distribution of GNPs. However, no effective imaging tool exists to determine the GNP distribution in vivo. Various approaches have been proposed to determine the concentration of GNPs and its distribution in a tumour and in other organs and tissues. X-ray fluorescence computed tomography (XFCT) is a promising imaging technique to do that. A new experimental device based on the XFCT technique allowing the in vivo control of GNP radiotherapy treatments is proposed. As a test of concept, experimental acquisitions and Monte Carlo simulations were performed to determine the performance that a XFCT detector has to fulfil.

INTRODUCTION

Contrast agents with high atomic number, such as gold nanoparticles (GNPs), have the potential to enhance the radiation dose locally in external radiotherapy with kV X-rays. In the kilovoltage region, gold has an attenuation coefficient two orders of magnitude greater than that of soft tissue. So, even a small amount of GNPs allows a dose enhancement inside the tumour volume by sparing surrounding normal tissue provided that the GNP concentration is lower there. The greatest ratio between the absorption coefficient of gold and that of soft tissue is found near 40 -50 keV. However, the use of such energies for treatment would be complicated by the high absorption in tissue. Therefore, X-rays of higher energies, near the gold K-edge (80.7 keV), are usually used. Several developments have already been performed with monoenergetic beams around the gold K-edge, using synchrotron radiation, but that technique is not useful in a clinical context. Instead, more traditional X-ray production methods are considered with some refinements in the filtration to optimise the photon energy distribution of the beam [START_REF] Ricketts | A quantitative x-ray detection system for gold nanoparticle tumour biomarkers[END_REF] .

One of the key tasks during preclinical studies of such a treatment method is to determine the biodistribution of GNPs injected into animals. Indeed, the K a fluorescent X-ray energies of gold (66.99 and 68.80 keV) are both able to escape a body when produced by polychromatic X-rays [START_REF] Popovtzer | Targeted gold nanoparticles enable molecular CT imaging of cancer[END_REF][START_REF] Lee | Functional CT: physiological models[END_REF] . However, it has been noted [START_REF] Cheong | X-ray fluorescence computed tomography (XFCT) imaging of gold nanoparticle-loaded objects using 110 kVp x-rays[END_REF][START_REF] Ricketts | An x-ray fluorescence imaging system for gold nanoparticle detection[END_REF] that the following three technical challenges are yet to be met satisfactory to perform in vivo quantification of GNPs:

-Goal 1: imaging of GNPs distributed within the tumour and critical organs, -Goal 2: quantification of the amount of GNPs present, -Goal 3: achieving system sensitivity capable of detecting GNP concentration typically found in tumours.

X-ray fluorescence computed tomography (XFCT) has shown the potential to overcome these challenges by providing GNP images as well as concentration information (4, 6 -8) . XFCT imaging of GNP-loaded phantoms has been performed using K X-ray fluorescence (K-XRF) in gold [START_REF] Cheong | X-ray fluorescence computed tomography (XFCT) imaging of gold nanoparticle-loaded objects using 110 kVp x-rays[END_REF] . However, relatively high GNP concentrations were tested by the system. Simulation studies comparing X-ray fluorescence CT and K-edge CT show that fluorescence imaging can provide a better sensitivity for very low tracer concentration [START_REF] Bazalova | Investigation of x-ray fluorescence computed tomography (XFCT) and K-edge imaging[END_REF] . The development of a new experimental device based on the XFCT technique and allowing in vivo control of GNP radiotherapy treatments is now proposed. As a test of concept, in a first step the X-ray tube and the detector used at the authors' institute for non-destructive testing were used. In parallel, Monte Carlo (MC) simulations performed with PENELOPE were used; on one hand, as an optimisation tool to enhance the experimental set-up, and on the other hand, as a planning system to estimate the GNP concentration needed to treat tumours.

MATERIALS AND METHODS

Experimental set-up

The proposed system (Figure 1) was configured according to the experimental set-up reported by MacMahon et al. [START_REF] Mcmahon | Radiotherapy in the presence of contrast agents: a general figure of merit and its application to gold nanoparticles[END_REF] and consisted of a polychromatic X-ray source (125 kVp, 300 mA), a 6 mm diameter collimator, a spectrometer system (XR 100T-CdTe, Amptek, Inc., USA), a multichannel analyser (Novelec SM-1024) and a phantom.

The XR 100T-CdTe was a thermoelectrically cooled 5 `Â 5 `Â 1 mm 3 CdTe diode X-ray detector and preamplifier. Its energy resolution was 1.5 keV (FWHM) at 150 keV and 0.6 keV at 59.54 keV ( 241 Am source). It had high detection efficiency up to 100 keV, which made it ideal for XRF applications. Its calibration was performed with a metallic sheet of pure gold (99 %).

The phantom was composed of two PMMA parallelepipeds of 5 cm thickness mimicking the healthy tissue and a cylindrical 'tumour' volume of 5 cm diameter, positioned between the two PMMA parallelepipeds and in the intersection of the incident and detected beams. The tumour volume was made of PMMA. It had a small hole into which a sheet of pure gold of 1 mm height and 25 mm diameter was inserted. The dimensions of the sheet were chosen to obtain a concentration of gold in the tumour volume equivalent to 10 mg Au g 21 of tumour.

Several strategies were applied to effectively minimise the unwanted background in the registered pulseheight distributions: (i) two 0.6-mm-thick copper filters were placed at the exit of the collimated beam to reduce the unnecessary X-rays with energies below gold K-edge energy (80.7 keV), (ii) the spectrometer was positioned at 1358 to the incident X-ray beam to minimise the impact of the tungsten peaks and of the Compton scattered X-rays in the region of the Au peaks, and (iii) a brass collimator with an opening of 6 mm diameter was placed to cover the sensitive element of the spectrometer in order to receive the directional beam while improving the signal-to-noise ratio.

A first measurement was acquired with gold and without any PMMA phantom to detect the gold peaks and to calibrate the spectrometer's energy scale. A second measurement was performed with the PMMA inserts and the tumour volume without gold. A third measurement was performed with the PMMA inserts and the tumour cylindrical volume with gold. The three measurements lasted 10 min at 300 mA.

MC simulations

Simulations were performed with the 2006 version of PENELOPE [START_REF] Salvat | A code system for Monte Carlo simulation of electron and photon transport[END_REF] . An RQR10 spectrum [START_REF] Oliveira | Uniformity and field size of filtered x-ray beams[END_REF][START_REF] Oliveira | On implementing reference radiations for calibrating and testing solid state detectors for dose measurements in radiology[END_REF] was simulated with a 0.958 aperture. The simulated geometry corresponded to the experimental set-up. A brass collimator with a 6 mm diameter opening was placed on the beam axis, and all the photons emerging it were registered to build an energy spectrum. The same brass collimator was also positioned at 458 from the beam axis, and the photons emerging it were also registered.

A second simulation was performed with the same parameters, except that the tumour cylinder was replaced by water.

RESULTS

Experimental results

Figure 2 shows the spectra obtained with and without any gold sheet in the tumour volume. Both spectra present a Compton scattering front between 62 and 66 keV. The K b peaks of tungsten are visible at 67.2 and 69.1 keV for the 'PMMA spectrum' and the Au K a peaks can be detected at 67.0 and 68.8 keV for the 'Gold þ PMMA spectrum'. The difference between both spectra is plotted in the second panel of Figure 2. Au peaks are hardly detected due to the high background level.

Simulated results

Spectra simulated with and without gold in the tumour volume are reported in Figure 3. The difference between the two simulated spectra is shown in the second panel of Figure 3. It shows a very good discrimination of the Au K a2 peak, due to the better energy resolution used in the MC simulation and to the absence of background induced by the detector. However, the Au K a1 peak is hardly detectable, probably being inferior due to the MC uncertainty.

DISCUSSION

This pilot study has been performed for one experimental and one simulated configuration. The results show that gold peaks are not easily detected because their contribution is mixed with tungsten K a peaks, phantom Compton scattering and interactions in the detector. It has been shown [START_REF] Mcmahon | Radiotherapy in the presence of contrast agents: a general figure of merit and its application to gold nanoparticles[END_REF] that an optimisation of the beam filtration can be performed to optimise the spectral shape of the excited source and increase the sensitivity. Moreover, the detector position has also to be tested in order to find the most suitable detection set-up. Scattered X-rays can be minimised by using a backscatter detector orientation [START_REF] Ahmad | Order of magnitude sensitivity increase in x-ray fluorescence computed tomography (XFCT) imaging with an optimized spectro-spatial detector configuration: theory and simulation[END_REF] . Further simulations will be performed to determine the experimental set-up that is most suitable for application in the present study. To improve the simulation, the detector model will be integrated into the simulation to analyse its impact on the Au peaks detection. Moreover, in a first step, the authors' laboratory detector usually used for applications in non-destructive testing was used to check the feasibility of the method. However, the energy resolution of this one is too bad (1.5 keV at 150 keV) to allow a good discrimination of the gold peaks in the Compton background. A better energy resolution is needed to allow a good discrimination of the Au peaks.

In the future, the authors plan to test the performance of a detector initially developed in their institute for Astrophysics [START_REF] Meuris | a CdTe based spectrometer for bright solar event observations in hard X-rays[END_REF][START_REF] Dubos | a CdTe-based spectro-imager development for nuclear applications[END_REF] . This is an imaging spectrometer (CALISTE HD), which uses a 1-mm-thick Schottky CdTe detector divided into 256 pixels disposed in a 16 `Â 16 matrix with a pixel pitch of 625 mm. It has been tested with a 241 Am source, and excellent spectroscopic performance was obtained, in particular between 10 and 100 keV (0.56 and 0.67 keV FWHM at 13.9 and 59.5 keV, equivalent to 4 and 1.1 % energy resolution, respectively). A better energy resolution ( 0.52 keV at 59.5 keV) can be attained for the present application after a finer calibration step and a reduction of the electronical noise.

CONCLUSION

Measurements performed during this study prove the feasibility to detect Au peaks in GNP radiotherapy treatments. The MC simulations show that a detector with a better energy resolution than the one used during this study is needed to get a good discrimination of the peaks. Other parameters of the detector have also to be taken into account, such as the calibration process, the scattering interactions into the elements of the detector, the pile-up effect and the dead-time to minimise the background visible in Figure 2. A more detailed MC simulation taking into account the geometry of the detector could provide useful information to characterise the ideal detector. Experimental acquisitions will thus be conducted with a better energy resolution detector developed initially for astrophysical applications. In parallel, MC simulations will be performed to optimise the experimental set-up and quantify the performance of the detector.

Figure 1 .

 1 Figure 1. Experimental set-up: the tumour volume is placed on the PMMA cylinder and between the two PMMA inserts and the incident beam is impinging perpendicularly on them. The spectrometer is placed at 1358 from the incident beam.

Figure 2 .

 2 Figure 2. In the first panel, experimental spectra obtained with (a) and without (b) gold inserted into the small cylindrical phantom are plotted. In the second panel, the difference between the two spectra is plotted. The dashed and dotted lines indicate the K a X-ray energies of gold (66.99 and 68.80 keV) and the dotted lines the K b X-ray energies of tungsten (67.244 and 69.1 keV).

Figure 3 .

 3 Figure 3. Simulated spectra obtained with (a) and without gold (b) inserted into the small cylindrical phantom. In the second panel, the difference between the two spectra is plotted. The dashed and dotted lines indicate the K a fluorescent X-ray energies of gold (66.99 and 68.80 keV).