Cyrille Balland
email: cyrille.balland@ineris.fr

Vincent Renaud
email: vincent.renaud@ineris.fr

High-resolution velocity field imaging around a borehole: Excavation Damaged Zone characterization

The excavation of a deep underground structure induces a stress field redistribution that creates an Excavation Damaged Zone (EDZ). Study of EDZ physical properties is of prime importance in the framework of nuclear waste geological storage, where EDZ constitutes a preferential pathway for storage materials towards the biosphere. Use of ultrasonic wave propagation inside the rock media is of great help as a near field non-intrusive technique for this. The significant research activity in the laboratory offers a large number of boreholes. These boreholes form linear excavations with a perfectly round section. It is for this reason that a specific probe and inversion method have been developed in order to reconstruct the image of the velocity field all around the environment of the boreholes. In order to recover all data needed with high accuracy and minimum effort, the probe is fully monitored for automatic data acquisition. This numerical inversion method has been preliminarily applied to a multi layer test bed in order to validate and optimise the experimental and numerical procedures. Based on those results, an in situ experiment has been performed in the Meuse/Haute-Marne Underground Research Laboratory (URL) in France for the safety assessment of nuclear waste storage in the Callovo-Oxfordian argillaceous layer.

INTRODUCTION

Rock damage introduced by Kachanov in 1958 corresponds to the development of microscopic cracks in the whole of a material. In the case of underground storage, rock damage will affect the rock capacity to confine radionuclides. The principal factors involved in initiation and extension of damage result from the excavation method, the shape of the opening, the stress field and the mechanical behaviour of the rock mass. The stress field and the geometry are generally known. However, the in situ rheology of the rock is much more difficult to determine. Its mechanical characteristics are often variable, particularly in the case of sedimentary rocks that have significant heterogeneities and high anisotropy depending on the settlement. In addition to the natural variability of the rock, there are also uncertainties linked to the experimental methods used to characterise the geomechanical properties of the rock. It is consequently difficult to set up a constitutive model and even more difficult to predict the level of damage that an excavation in a given rock medium may generate.

In order to study EDZ, ultrasonic monitoring and auscultation have proved to be efficient [START_REF] Maxwell | Seismic imaging of rock mass responses to excavation[END_REF]. Deformations and microcracks markedly disrupted the propagation of ultrasonic waves [START_REF] Nur | Stress-induced velocity anisotropy in rocks: an experimental study[END_REF]. Simple ultrasonic devices are able to detect the initiation of microcracks at a very early stage. When these devices process a larger amount of information, it is then possible to reconstruct an image of the medium by ultrasonic tomography. Tomographic images have already been obtained in underground structures in order to investigate the damaged zone [START_REF] Maxwell | Seismic imaging of rock mass responses to excavation[END_REF][START_REF] Forney | Caractérisation par méthodes ultrasoniques de la zone endommagée induite par le creusement d'un tunnel en milieu argileux : Cas d'étude au tunnel du Mont Terri[END_REF][START_REF] Damaj | Velocity survey of an excavation damaged zone: influence of excavation and reloading, Rock Physics and Geomechanics in the Study of Reservoirs and Repositories[END_REF]. In order to have a good resolution, however, it is necessary to deploy significant resources with a large number of sensors and boreholes [START_REF] Balland | Ultrasonic velocity survey in Callovo-Oxfordien argillaceous rock during shaft excavation[END_REF].

Instead of deploying considerable efforts, an alternative described in this paper consists in studying the damaged zone induced around a borehole drilled for survey purposes. This would involve the installation of several boreholes around big-sized excavation and reconstruction of the complete image of the damage around each borehole. With the exception of the scale effect (from the borehole to the structure), the results of this procedure are expected to more accurately estimate the damage, and more generally the behaviour of the structure.

Previously, [START_REF] Hornby | Tomographic reconstruction of near-boreholes slowness using refracted borehole sonic arrivals[END_REF] has demonstrated that it is possible to reconstruct the field in a longitudinal section along a borehole. Winckler (1996) has highlighted azimuthal velocity variations due to the stress anisotropy and has developed an ultrasonic imaging method of the wall with D' Angelo (2006). In our study, the developed method consisted in reconstructing the image of the velocity field around a borehole in a radial plane. This method has recently been developed and validated within the framework of the European NF-Pro project. It has been applied for the first time in the Meuse/Haute-Marne Underground Research Laboratory. It is based on an innovative ultrasonic probe capable of automatically performing a large number of measurements in order to build a "high resolution" radial image of the velocity field around a borehole. This paper is composed of the 3 principal stages in the order of the study. The first stage is purely numerical with the development of an inversion method of the velocity field around a borehole. The second stage was to manufacture a test bed, choose the necessary instrumentation and test the inversion method. The last, but not least stage was to design and manufacture a probe tool and apply it in the site.

MEASUREMENT METHOD PRINCIPLE

The disturbed velocity field around a borehole has a mainly radial gradient. The lowest velocity is generally close to the wall and increases progressively with depth (Figure 1). A stress wave emitted and received from the surface of the borehole will tend to probe the zones of the rock with the highest velocities, thereby respecting Fermat's principle. Then with just a pair of emitter and receiver sensors whose relative distance extends linearly, it is thereby possible to probe the rock deeper and deeper. Once radial coverage is achieved, a slight rotation of the probe around the borehole axis and moving it along the borehole, the whole volume of the surrounding rock can be monitored, providing reconstructed 3D information on the P-wave propagation velocities.

INVERSION METHOD

The inversion process of ultrasonic data to velocity field images is based on a linearized iterative reconstruction technique significantly dependant on its initial solution. Although rule-of-thumb considerations are commonly used and tested to fix this initial solution, a different refined approach was undertaken by the authors.

It is based first on numerical work on the calculated stress field surrounding a circular shaped excavation taking the calibrated complex rheology of the argillites into account [START_REF] Souley | Constitutive models for rock mass numerical implementation, verification and validation[END_REF] and, secondly, on empirical relationships established between stresses and velocities of ultrasonic waves from tests run on rock samples in laboratory conditions [START_REF] Pettitt | Development of the Tools and Interpretation Techniques for Ultrasonic Surveys to Monitor the Rock Barrier around Radioactive Waste Packages in Geological Repositories[END_REF]. It appeared worthwhile to obtain a first approximate order of the initial solution by changing the stress field by the corresponding P-wave velocity field around the circular borehole.

More precisely, the velocity field modelled (Figure 2d) is based on empirical relations [START_REF] Renaud | Numerical simulation and development of borehole ultrasonic imaging data inversion[END_REF] linked to the deviatoric stress (Figure 2b), the minor stress (Figure 2c), and the level of damage of the rock: elastic, damaged or fractured (Figure 2a). This calculated velocity field model is assumed to be close to that expected in situ. However, it should be noted that the empirical relationship between velocity and stress variations is based on very little experimental data and under load conditions not exactly representative of the stress state expected around the excavation.

The second step of the inversion method consisted in defining, as accurately as possible, the geometry of the rays' path in the planes along the directions of the major and minor horizontal stresses. To do this, modifications were made to improve a ray tracing algorithm originally developed for analysing earthquakes [START_REF] Um | A fast algorithm for two-point seismic ray tracing[END_REF]. Offset coefficients are optimised so that each of the rays probes the zones of the velocity field that minimise the travel time. This method has proved to be effective in most cases, except when the velocity gradient is too high locally (>100 m/s per mm). It is therefore necessary to check the validity of the method for high velocity gradients including a safety margin (>50 m/s per mm).

The inverse problem defined by the velocity field and ray path can be solved by a tomographic inversion, also known as global matrix inversion. Several other inversion methods were explored with different devices and different ray coverage. A Monte-Carlo method and another method based on a step-by-step reconstruction have demonstrated good precision [START_REF] Renaud | Numerical simulation and development of borehole ultrasonic imaging data inversion[END_REF]. However, the global matrix inversion proved to be the fastest method with a sufficient resolution close to the uncertainties of P-wave picking and the sensors' position. The quality of the inversion depends mainly on the quality of the experimental data. It is also possible to improve the overall result by optimising the damping coefficient [START_REF] Bois | Essai de détermination automatique des vitesses sismique par mesure entre puits[END_REF]. In our case, this involves minimising the variance on the travel times at the same time as the variance on the slowness residuals. The coefficient is optimal when it reaches the knee point of the curve of the variance of path time residuals as a function of the variance of slowness residuals. All of these new developments, as well as the principal ultrasonic tomography functionalities, have been developed in the TOMSIS © software.

NUMERICAL VALIDATION

The inversion method previously described was tested for the two profiles from velocity modelling around the borehole. These two sets of synthetic data also made it possible to test the two possible ultrasonic configurations from the wall of a borehole with configurations 1 and 2 (Figure 3) with rays centred respectively along a single median and all of the possible rays according to a constant displacement of the sensors.

The inversion grid for 0.08 m in depth is radially composed of a cell every 3 mm with 50 different sensor positions. There are therefore 25 rays in configuration 1 and 1225 in configuration 2. The travel times are calculated according to the optimised rays' paths (Um and Thumber, 1987 algorithm improved by [START_REF] Renaud | Numerical simulation and development of borehole ultrasonic imaging data inversion[END_REF] in the two velocity fields (major and minor stress) and with the two configurations. To these travel times is added a random noise of 0.5 μs corresponding to the approximate error on the picking of the first arrivals.

The data collected is then inversed following an initial velocity field with constant gradient. After the first inversion, the rays are replotted along the calculated velocity field. The optimal damping coefficient is recalculated for the second inversion. At the end of the third inversion, the variance of the residuals no longer changes significantly. Further inversions do not significantly improve the results. These results are shown in Figure 3, where the velocity difference between the initial velocity field and the inversed velocity field does not exceed 120 m/s except for the zone of the highest gradient on the profile in the direction of the minor stress. This confirms the limits of the method related to the amplitude of the velocity gradients. The variations of the velocities as a function of the initial stress orientation and the depth are clearly highlighted, and validate the overall method. Several factors, such as the number of symmetrical rays or the initial velocity field, could improve the inversion results.

Configuration 2, which has a much higher density of rays than configuration 1, does not significantly improve the resolution of the velocity profile, mainly because the velocity does not vary longitudinally to the borehole. The additional difficulty of implementation is then not justified here. It could be interesting to use configuration 2 to detect longitudinal heterogeneities as strong sources of errors. The method is therefore focused on configuration 1 because configuration 2 is much more technically complex to carry out and the argillites are quite homogeneous in the probed rock volume.

PHYSICAL MODEL

The synthetic study presented above shows the acceptable resolution and stability of the inversion of travel times collected to image the velocity field and the EDZ. To refine the validation and testing of the described approach, a physical model was built up to experiment data acquisition and processing as close as possible to the expected field conditions.

DESIGN

In order to reconstruct a medium with the velocity profile around the borehole, a precision of about a few m/s on the P-wave velocity is necessary. To do this, the physical model was manufactured based on a mixture of cement and clay (kaolinite) so as to have a final material with mechanical characteristics quite similar to those of argillites. The cement grants the mixture high mechanical strength and low attenuation, whereas the kaolinite fraction introduces a gradual reduction of both the mechanical strength and the P-wave velocity while increasing the attenuation. Nevertheless, acoustic measurements carried out on a variable mixture of kaolinite and cement show that the P-wave velocity of cement is greater than that expected in the target rocks around 5200 m/s for 3500 m/s in argillites, introducing an important P-wave component (Figure 4). In any case, only the relative variations of the Pwaves' velocity between the layers estimated during the modelling were therefore respected. In practice, this meant shifting the velocities of the numerical model while at the same time remaining positioned in a range of velocities appropriate for the cement/kaolinite mixture.

The multilayer test bed was ultimately made up of 17 layers from 3 to 110 mm thick, representing a total thickness of 203 mm (Figure 5). The number and thinness of the layers were directly proportional to the velocity changes. The thinner layers are close to the surface, where the velocity gradient is steeper.

LABORATORY MEASUREMENTS

The measurements were therefore taken from the surface of the test bed corresponding to the final layer poured on, with the lowest velocity. The measurement device was composed of 2 remote sensors, with a distance (± 0.1 mm). The first sensor fixed to the test bed acted as an emitter with a pulsed signal centred on a frequency of 65 kHz, while the second sensor was displaced on the test bed with a pitch of 1 to 2 mm (Figure 5). The emitter and the receiver sensors were connected via an amplifier to an ultrasonic acquisition unit with a sampling frequency of 10 MHz.

The analysis of the propagation of the waves in the test bed was possible up to an inter-sensor distance of 0.365 m. Beyond this value, the wave attenuation was too high: it led to a signal to noise ratio (<3) that did not enable a manual travel time picking to be made or crosscorrelation of the signals between each other. The apparent velocity measured on the profile (Figure 6: P-wave apparent velocity curve as a function of the inter-sensor distances for the multilayer test bed. Figure 6) varied from 3350 m/s for a distance between the two sensors of 0.03 m to 4900 m/s for a maximum inter-sensor distance of 0.365 m. The increase in the apparent velocity was globally monotonous with an asymptotic evolution for the longest rays (>500 m/s). The inversion of this data set was therefore possible according to a bending ray and global tomography approach. The velocity model for the inversion of the profile is composed of 19 layers of 0.003 m each, which form a grid 0.057 m deep. The thickness of this model chosen for the inversion corresponds to the maximum investigation depth determined by a sensitive study. The optimum number of layers constituting the model was set in relation to the number of rays. It is necessary for each cell to be crossed by a sufficient number of rays with sufficient azimuthal coverage so that the velocity affected is properly constrained. From a numerical point of view, the inversion matrix becomes a sparse matrix, and therefore more difficult to inverse with a slower convergence. Beyond 19 layers, the resolution of the final velocity field was no longer enhanced.

The results of the inversion are shown in Figure 7 where the "apparent" velocity is plotted as a function of the depth. The investigation depth of the experimental data corresponds to the maximum depth of the tomographic grid (red curve and triangles). The distances corresponding to the intermediate cells are estimated by simple linear interpolation. The investigation depths are not realistic at that stage, but they are of the same order of magnitude as those of the inversion making it possible to visualise the profile of inversed velocities (green curves and diamonds) at the same time as that of the apparent velocities measured on the test bed. On the same graph, the actual velocity profile of the test bed is shown (black curve and crosses). The inversed velocity profile stops a little before 0.04 m, which corresponds to the depth of the ray that explores the deepest layer. The inversed velocity profile is distinguished from the test bed profile by two main aspects: the velocity is lower by 600 m/s in the zone near the surface (<7.5 mm), and the velocity of the inversed profile is higher by 200 m/s between 7.5 mm and 34 mm, even though the shape of the two curves is very similar in this zone.

The asymptotic and homogeneous shape of the apparent velocity variation profile as a function of the distance between the sensors shows that the test bed (Figure 6) constitutes a simulation representative of the expected velocity field around a borehole. The profile does not show rectilinear zones that are typical of refracted waves. This shows that the number of layers and the velocity variations between the layers are sufficient to sample the evolution of the velocity gradient. The first wave recorded may therefore reasonably be considered a pressure wave always transmitted in velocity fields with vertical and positive gradient. The test bed was therefore able to fulfil its first objective of setting and validating the inversion method for the apparent velocities. However, the results of the inversion differ with the initial velocity model expected in the test bed, particularly with velocities that were much lower than expected close to the test bed surface (-16% or -700 m/s). This difference stabilised around 200 m/s deeper than 0.01 m. The residual travel times of the tomographic inversion that lead to several tens of metres per second at most on the velocities cannot explain these differences. These velocity differences are therefore related to the manufacture of the test bed. Two explanations may be put forward:

(1) During the drying of the cement/kaolinite mixtures, the thinnest layers (3 mm) and those with the highest kaolinite contents showed visible cracks on their surfaces in contact with air. These cracks could explain the much lower velocities than expected in the superficial layers. This phenomenon is more important for the thinnest layers and would explain why it diminishes with depth, where the thickness of the layers increases.

(2) To explain the velocity difference on the rest of the profiles (beyond 7.5 mm depth), one may invoke the effect of a transverse anisotropy induced by the tabular and multilayer structure of the test bed. Anisotropy is notable with this type of structure, in sedimentary rocks, particularly argillites in eastern France. If the velocity difference is entirely due to a transverse anisotropy and considering the base velocity as the theoretical and measured profile plateaus (respectively 5100 and 5300 m/s), then the anisotropy factor that could explain this difference would be around 4%, in agreement with the level of anisotropy that can be found in this type of material [START_REF] Darot | Acoustic wave velocity and permeability, evolution during pressure cycles on a thermally cracked granite[END_REF]. The anisotropy effect is probably added to that of cracks for the superficial layers, which would reduce the drop in velocity due to the fissuring to 500 m/s. The error of 4 to 16% between the experimental data and the p model does not enable the developed inversion method to be completely validated. To do this, it would be necessary to quantify the impact of the surface fissuring, which is quite difficult given the thickness of the layers and the transverse anisotropy created by the layers. We did, however, continue this approach on-site, keeping this possible error in mind.

ULTRASONIC PROBE DESIGN

The results from the velocity field simulation [START_REF] Renaud | Numerical simulation and development of borehole ultrasonic imaging data inversion[END_REF]) and the physical model experiment were used directly to qualify the feasibility of an application of the method. In particular, those results helped to determine the ratio between the investigation depth and the distance separating the two sensors. This defined the sufficient precision in terms of travel time and location of the sensors on this basis. An ultrasonic probe was specifically developed for in situ ultrasonic velocity measurements (Figure 8). It is made up of 3 principal elements. The first (A) is fixed and comprises amplification components. A device for rotating a second element of the probe (B) is also located at this point. This part itself is equipped with a device for the translation of a module of two ultrasonic sensors (C). This third part has a pneumatic system applying a mechanical coupling of the sensors against the wall of the borehole, thereby allowing the probe to be used in dry boreholes. The ultrasonic sensors operate on a frequency band from 50 to 200 kHz. Data acquisition is entirely automated, resulting in a high precision in the positioning of the sensors both in rotation and in translation with ± 0.5° and ± 0.2 mm, respectively. The distance between the sensors can vary between 2 and 600 mm. This ultrasonic probe and the associated auscultation method were used in a vertical borehole drilled for survey purposes in one of the galleries of the Meuse-Haute/Marne URL (Figure 9). This gallery is oriented in the direction of the minor horizontal stress (N65). The vertical stress at this depth was estimated at 12.7 MPa (12.4 MPa in the direction of the minor horizontal stress). The anisotropy is 1.3 between the minor horizontal stress and the major horizontal stress [START_REF] Delay | Making of underground scientific experimental programme at the Meuse/Haute-Marne underground research laboratory[END_REF]. The measurement device was placed 5.10 m below the gallery level, where the influence of the gallery was sufficiently low. The homogeneous facies at this depth (-495.1 m) corresponds to more or less silty argillites with a Young's modulus of about 4 GPa and a uniaxial compressive strength of nearly 12 MPa. The experimental measurement plane was composed of 10,728 measurement points over the full horizon (360°) with a measurement pitch of 5° and between 37 mm and 337 mm for the inter-sensor distance with a pitch of 2 mm.

Initially according to straight rays with an apparent velocity, the results will be now corrected following the inversion of arrival times to go up at a "realistic" velocity following the bedded rays. The non-corrected results of the different measurements on the multilayer test bed and in the borehole are given in the form of apparent velocity curves as a function of the inter-sensor distance. The first arrival of the incident wave corresponds to that of the P-wave. The first picking is carried out manually on the waveform in which the signal to noise ratio is the highest. The other waveforms are then compared to this reference signal by cross-correlation by threshold overrun. The set of data was then inversed by ultrasonic tomography according to the method described above.

BOREHOLE MEASUREMENTS

An acquisition methodology similar to that of the physical model was applied with the ultrasonic probe in the borehole. In the case of the test bed, it involved only one plane perpendicular to the layers. For the borehole, it involved several radial planes passing through the axis of the borehole. The problem can be considered as similar due to a velocity gradient expected with a single radial component, which varies with the depth and the azimuth. To cover the horizon of the borehole as well as possible, 72 measurement semi planes were arranged with a 5° step (Figure 10). In each semi plane, 150 P-wave propagation time measurements were carried out. Another part of the data was ruled out, due to the overly low signal to noise ratio, which did not make it possible to pick the first arrival particularly for the complete N195° profile and as well as the part of data with significant propagation lengths. Finally, out of the 10,728 measurements, 8,390 measurements were retained for the analysis, i.e. nearly 78%.

All of the raw data are presented in Figure 11 and Figure 12. Each grey cross gives the "apparent" velocity calculated on the straight ray corresponding to one measurement configuration. The points of measurements are present according to the inter-sensor distance. The four profiles corresponding to the principal directions of the stresses are presented differently (blue rhombuses and red triangles in the direction of the horizontal minor stress, pink square and black circle in the direction of the major horizontal stress). In Figure 12, the points of measurements are presented according to the horizon.

The "uncorrected curves" show that the "apparent" velocity (Figure 11 and Figure 12) is minimal with 1630 m/s for an inter-sensor distance of 0.041 m and an orientation N60°. Its maximum is 2705 m/s for a distance of 0.313 m and 340°. The evolution of the velocities for different orientations was greater (400 m/s) for the shortest inter-sensor distances (0.037 m) than at the end of the profiles (200 m/s) when the inter-sensor distance was highest (0.337 cm). The profiles showed an asymptotic increase of the same type as those of the numerical and physical models with, however, a break of slope for inter-sensor distances around 0.013 m for profiles oriented north-east. Velocity variations as a function of the angle of the measurement plane were more significant with the shortest inter-sensor distances than with the farthest away sensors. In Figure 12, for example, the velocity curve with a spacing between the sensors of 0.037 m showed velocity variations around the horizon of more than 300 m/s, whereas the curve with 0.237 m shows variations of around 100 m/s. For the shortest inter-sensor distances, the variations were abrupt and connected the two "velocity plateaus" by less than 20°; only several measurement points made the transition between these azimuthal spaces that form the plateaus where the velocity varies little. These azimuthal velocity profiles showed a very pronounced symmetry of revolution for all inter-sensor distances. In Figure 11, the profiles N65° and N245° corresponding to the direction of the minor horizontal stress are relatively well superimposed; the same holds true for the profiles oriented N155° and N335°, which correspond to the direction of the major horizontal stress. The apparent velocity difference between these two groups of profiles is significant when the sensors are very close, and then decreases rapidly when the sensors are farther apart.

All of these profiles may be represented on a velocity map of the investigated section (Figure 13). This figure presents the interpolating image of all "apparent" velocities (also presented in Figure 11 and Figure 12) function of the azimuth and the geometrical radial distances of the measurements. This image assumes that the velocity variations are mainly radial and azimuthal and that they are negligible in the axis of the borehole, over the maximum distance separating the two sensors. Acoustic logs carried out in this zone [START_REF] Vinsot | Référentiel Géologique du site de Meuse/Haute-Marne[END_REF] have shown that the velocity has a very low vertical variability of about 100 m/s per m. The variation between the two extreme distances is therefore close to 30 m/s, which may be disregarded with respect to the measured radial variations.

Moreover, due to the anisotropy, the "true" velocity (unlike the "apparent" velocity) measured along the ray path results from the combination of velocities in the different directions that the ray takes, radial near the sensors and longitudinal midway. It is therefore assumed that this combination has little influence on the final measurement. The relatively low anisotropic velocity of argillites, naturally 10%, decreasing with damage [START_REF] Pettitt | Development of the Tools and Interpretation Techniques for Ultrasonic Surveys to Monitor the Rock Barrier around Radioactive Waste Packages in Geological Repositories[END_REF], and the shape of the rays (Figure 14) gives a weak anisotropic effect on the travel times. This effect has therefore not yet been taken into account in this study.

In Figure 13, the localisation of the measurement points is fixed arbitrarily regarding the depth of the inversion grid while maintaining the same proportions between the different pairs of sensors, in the same way as for the apparent velocity profile on the test bed (Figure 6). In other words, the depth given in this graph is only dependent on the distance between the pair of sensors. This representation clearly shows symmetry of revolution for the apparent velocities. It also appears that the velocity profiles are more heterogeneous in the directions close to the minor horizontal stress than in the direction of the major horizontal stress, probably caused by the stronger damage in this direction.

The previous numerical grid for 0.057 m in depth composed of 19 layers is therefore still valid. The investigation depth reached by the rays was, however, variable as a function of the orientation of the planes due to the fact that the calculated velocity field differs significantly between the various planes but always in the domain of 0.057 m. Figure 13 shows the profiles of apparent velocities in all the monitored directions and in the directions of the minor and major horizontal stresses, as well as the four corresponding inversed profiles. For these profiles, the investigation depth was similar, corresponding to the depth of the centre of the last tomographic cell crossed by the rays (Figure 14). Since the velocity gradient is positive (with increasing depth), these rays are logically those that have the greatest inter-sensor distances.

The inversed velocities varied from 1642 m/s at 0.0015 m from the wall for an orientation of N65° to 3029 m/s at 0.0255 m from the wall for an orientation of N155° (Figure 15). As for the apparent velocities, the scattering of the inversed velocities was greater inside the immediate wallside (>400 m/s) and was reduced in greater depth (<200 m/s). This scattering was, however, greater than that of the apparent velocity profiles. The shape of the curves showed a more complex variation of the velocity gradients with moderate gradient near the wall, to the steepest just after 0.01 m depth, then less and less until becoming negative in the zone furthest from the wall. This inversion occurred more or less soon depending on the profile. The four profiles taken in the directions of the principal horizontal stresses confirm the symmetry already observed in the apparent velocity profiles. In fact, profiles of the same orientation do not show any significant difference except for the two profiles N155 and N335, where the velocity difference reached nearly 80 m/s at a depth of 0.045 m. Between the profiles with different principal stress directions, the difference in velocity was maximal (from 440 to 520 m/s as a function of the profiles chosen) near the wall, becoming negligible at 0.0105 m where the curves become tangential. Beyond this, the difference again increased, reaching a maximum of 190 m/s at 0.0195 m in depth, then diminished slightly at the end of the profile. As well as the profiles of the measured apparent velocity, the inversed velocity profiles are interpolated between them to constitute a 2D cross section map of the velocity field of the section of the studied borehole (Figure 16), this time with realistic velocities and depths. This image shows a crown of low velocity (<2500 m/s) of about 0.02 m thickness all around the borehole, with a zone of very low velocity (<2000 m/s) of about 0.01 m in depth localised at the wall in a first zone between N10° and N95°, and in the opposite side between N200° and N295°. In these same angular zones, the isovalue line 2500 m/s was also more perturbed than elsewhere. With increasing depth, the velocity was globally less than the value in the other complementary angular zones.

DISCUSSION

Borehole measurements showed the same type of profiles for the physical model. The natural variations in the P-wave velocity at this depth (-495.1 m) are about 100 m/s per m over the full height (0.337 m) of the ultrasonic device [START_REF] Vinsot | Référentiel Géologique du site de Meuse/Haute-Marne[END_REF]. The maximum investigation depth is close to 0.05 m for the maximum inter-sensor distance. The anisotropy of the P-wave velocity in the medium is transversal and does not exceed 10% [START_REF] Balland | Ultrasonic velocity survey in Callovo-Oxfordien argillaceous rock during shaft excavation[END_REF]. The influence of the anisotropy on the rays will therefore be 1.5% at most. Before excavation of the borehole, and given the expected perturbations, the medium may be considered acoustically homogeneous in the monitored volume and isotropic in the axis of the device. The majority of the velocity variations measured therefore results from the perturbations induced during the excavation of the borehole. The variation in velocity as a function of the investigation depth reaches a maximum of 60% to 70%, respectively, for the apparent velocity and the inversed velocity. In other words, there is an increase from 60% to 70% between velocity at the wall and the deepest point observed by the device. This variation is only 30% to 35% in the direction of the major horizontal stress. These velocity variations correspond to: (1) the variation of the effective stress field and the elastic characteristics of rocks in general [START_REF] Nur | Stress-induced velocity anisotropy in rocks: an experimental study[END_REF], and in particular for argillites as reported by [START_REF] Pettitt | Development of the Tools and Interpretation Techniques for Ultrasonic Surveys to Monitor the Rock Barrier around Radioactive Waste Packages in Geological Repositories[END_REF] (2) a coalescence of microcracks induced by the deviatoric stresses [START_REF] Lockner | Changes in seismic velocity and attenuation during deformation of granite[END_REF][START_REF] Hudson | Wave speeds and attenuation of elastic waves in material containing cracks[END_REF][START_REF] Schubnel | Velocity measurements and crack density determination during wet triaxial experiments on Oshima and Toki granites[END_REF]. These phenomena are also accompanied by perturbations of the saturation and the porosity of the rock in the near field [START_REF] Murphy | Effects of partial water saturation on attenuation in Massilon sandstone and Vycor porous glass[END_REF][START_REF] Krief | A petrophysical interpretation using velocities of P and S waves[END_REF] as well as less significant phenomena. At this stage of development of the method, the measurement of only the P-wave velocity does not enable us to separate the exact contribution of these different phenomena in the perturbations of the propagation of the waves.

The P-wave velocity, which decreases less in the direction of the major principal stress than in the direction of the minor stress, is close to what has been observed in the laboratory on sandstone specimens taken from a borehole (Winckler and D'Angelo,, 1996) where the Pwave increases in the direction of the major principal stress and decreases in the direction of the minor stress. The difference of the tendency in the major principal stress for the sandstones and argillites comes from the mechanical behaviour difference. The sandstone's mechanical behaviour is quite elastic with anisotropic damaging initiated close to the break point. Quite the reverse, the argillites' damaging occurs sooner, just after a short elastic phase, and is quasi isotropic. Indeed, the sandstones' loading could generate some local velocity increases, unlike the argillites which quickly sustain plastic strain, damaging and velocity decrease. The other notable difference is that the velocity variation is much more abrupt in argillites, particularly for the rays closest to the borehole wall. As has been seen previously in Figure 11, Figure 12, and Figure 13, the evolution of the velocity is quasi-discontinuous, which is contrary to the progressive variations of the stress field around a cylindrical excavation. This again can be explained by the high damaging or fracturing observable in argillites. The fracturing of the rock at the surface of the wall forms the well-known "breakouts" already observed in different underground laboratories [START_REF] Martin | Observations of brittle failure around a circular test tunnel[END_REF][START_REF] Thury | The Mont Terri rock laboratory, a new international research project in a Mesozoic shale formation[END_REF]. This phenomenon may also be correlated with the visual observation of the experimental borehole (Figure 17). Traces of core drill on the wall showed an out-ofroundness from 3 to 4 metres deep and a damaged wall state in the direction of the gallery (direction of the minor horizontal stress). The out-of-roundness and the "breakout" only appear at a certain depth, which may be explained by the influence of the gallery on the stress field. Near the gallery, the stresses applied on the borehole wall and the stress anisotropy are minimised.

The borehole measurements also show that the asymmetry of the velocity variations diminishes with the monitoring depth (Figure 16). The velocity field goes progressively to the homogeneous and natural velocity in the rock. At the same time, and still with depth, the velocity transition between the velocity zone in the direction of the minor stress and that in the direction of the major horizontal stress becomes less and less abrupt. The velocities measured at this depth correspond to the higher inter-sensor distances. Moreover, the measured velocities integrate a greater rock thickness. It is certainly for these reasons as well that the azimuthal variations in the velocity measured in the other experiments are less abrupt.

It is also interesting to compare the image of the inversed velocity field with the results of the numerical model around the borehole proposed by Renaud in 2006 (Figure 2). Firstly, the amplitude of the velocity variations is minimised in the model with variations from 2200 to 3200 m/s and 1700 to 3000 m/s for the image of the inversed velocity field. In the direction of the major stress, the minimum velocity is now only 2200 m/s (Figure 15 and Figure 16), that is to say, the same order of magnitude as the numerical model. This difference could be explained by the size of the grid cells in the model, which cannot change the phenomena detail at the immediate wall. Except for the zone very close to the wall, the zone of very high velocity gradient (from 2200 to 2900 m/s for the numerical model and from 1700 to 2900 m/s for the experimental results) coincides relatively well with a thickness of about 0.015 m (0.35 times the radius of the borehole). The numerical model results show that the shape of the velocity field resembles the field of the minor horizontal stress (which may be correlated with the hydrostatic field) rather than the distribution of the deviatoric stress or the damage. In those places where the confinement is considerable, the microcracks, able to develop, remain closed and do not influence the P-wave velocity significantly. However, in the unconfined zones, the damage, or even failure leads to a significant reduction of the velocity. The lowest velocity zones (<2000 m/s) in the experimental velocity field correspond in terms of amplitude and extend to the lowest numerical velocities (<2400 m/s) and the fractured zone of the numerical model.

Finally, the characterization of the velocity field is efficient to characterize the EDZ of the borehole as a complement to numerical modelling. However, the method has not yet been used to estimate the EDZ of a gallery. In this next stage, it will act for a large gallery, to take measurements at several drillings with the same orientation and in each facies crossed by the gallery. In theory, some is the law of behaviour (elastic or elastoplastic), the ratio between the radius of the excavation and the radius of damage is constant. It is then possible to extrapolate the amplitude and the extension of the borehole EDZ to that of the gallery itself. This extrapolation nevertheless requires few conditions: (1) the characteristics of the rock mechanics must be homogeneous on the entirety of the volume disturbed by the gallery on a survey scale like that of the gallery, (2) the constituent elements of the rock (grains, matrix, interface) must have a size smaller than the resolution of the ultrasonic method. Even if the argillites of the callovo-oxfordien are weakly heterogeneous, these two conditions are never totally reached. It is thus necessary for the states of stress at the wall of the gallery and the boreholes to be in the same behavioural field. Extrapolation would be impossible for a borehole wall under the damage threshold and a gallery wall over the damage threshold. It would ideally be necessary to apply the characterization method for various borehole radiuses and thus to check the linearity of the function between the borehole radius and the damage radius on the scales of the study field.

CONCLUSIONS

The aim of this study was to see whether it is possible, for EDZ characterisation, to quantify the velocity field around an in situ borehole. An ultrasonic method based on the reconstruction of tomographic images and an automated ultrasonic probe was developed. The probe was designed for the acquisition of velocity data at depth based on the distance between a pair of sensors and the circumference of the borehole through a rotation of the device. The automation of the tool assured both high precision in the localisation of the measurements and a high number of samples that are taken into account during reconstruction of the image of the velocity field. This method was applied to a multilayer test bed performed to validate the method and to dimension the probe used for the first time in the ANDRA URL at -490 m. The results show that the method is capable of reconstituting a velocity field both in terms of amplitude and localisation. The velocity amplitude was not, however, faithfully evaluated by the method, particularly for the test bed, unlike the velocity field contours, which are in good agreement with the velocity field obtained by both numerical modelling and laboratory tests. The anisotropy of the stresses and the extension of the EDZ could thus be finely characterised. This method could therefore contribute to inferring the orientation of the principal stresses and the extension of the different perturbations of the rock properties in the immediate vicinity of the borehole.

The method can be applied to complete imaging around a borehole by displacing the measurement device along the borehole. However, the significant number of measurement points chosen for this study imposes an overly long measurement time. Moreover, the method does yet not take into account the velocity anisotropy of the rock medium. In the case of measurements with a direction anywhere in space, it would be necessary to measure and correct the influence of the anisotropy, to be able to compare each measurement plane and to obtain the overall velocity field of the rock mass. In the case of a more heterogeneous medium, it would be interesting to inverse a single matrix comprising all of the tomographic panels to best assess the longitudinal variations of the velocity field.

The approach of correlating the P-wave velocity considering the evolution of the rock properties shows the interest of the approach, in particular for a relatively homogeneous rock medium. However, the method is only based on one wave propagation parameter. The calculation of equivalent dynamic modules together with the measurement of S-waves would enable a better comparison of experimental data with the numerical models, with a view to better model calibration. This method could also contribute to predicting the damage around the excavations based on the velocity field both in boreholes with different directions and in the different rock.

This method can also be directly applied when the objective is to know the damage of a borehole such as, for example, within the framework of forthcoming thermal experimentation in the URL. It is necessary in this case to know the EDZ of borehole in order to take it into account in measurements of the thermal diffusion process. This method could also be used in the probing and the security of drillings for access to gas stores (hydrocarbons or CO2). Figure 3: Comparison between the velocity profiles of initial synthetic models in the directions of minor and major horizontal stresses and results of inversions using two different tomography configurations [START_REF] Renaud | Numerical simulation and development of borehole ultrasonic imaging data inversion[END_REF]. config.

1 1 1 1 2 2 2 2 config.

config.

Figure 3: Comparison between the velocity profiles of initial synthetic models in the directions of minor and major horizontal stresses and results of inversions using two different tomography configurations [START_REF] Renaud | Numerical simulation and development of borehole ultrasonic imaging data inversion[END_REF].

Figure 2 :

 2 Figure 2: Numerical modelling results with: a) The level of rock damage; b) The deviatoric stress; c) The minor stress and d) The P-wave velocity.

Figure 4 :

 4 Figure 4: P-wave velocity as a function of the radial distance for the numerical model (blue triangles) and the manufactured multilayer test bed (red crosses).

Figure 5 :

 5 Figure 5: Superimposed layers of the test bed and measurement device on the surface.

Figure 6 :

 6 Figure 6: P-wave apparent velocity curve as a function of the inter-sensor distances for the multilayer test bed.

Figure 7 :

 7 Figure 7: Comparison between the apparent velocities, expected and resulting from the inversion.

Figure 8 :

 8 Figure 8: Schematic diagram of the UPTO ultrasonic probe developedS.

Figure 9 :

 9 Figure 9: Meuse / Haute-Marne URL and experimental borehole location (courtesy of Andra (on left); Installation of the borehole probe (on right).

Figure 10 :

 10 Figure10: Layout of measurements around the borehole (the radial lay out of the measurements is not realistic, it corresponds exclusively to the projected inter-sensor distances) with, in red, the discarded measurements.

Figure 11 :

 11 Figure11: Apparent velocity curves of P-waves as a function of the inter-sensor distances for all of the measurements and in the direction of the principal minor and major horizontal stresses.

Figure 12 :

 12 Figure 12: Apparent velocity curves of P-waves as a function of the azimuth for all of the measurements and for different inter-sensor distances.

Figure 13 :

 13 Figure13: Image of the apparent velocity field (the depths are not realistic and have been defined in relation to the size of the tomographic inversion grid).

Figure 14 :

 14 Figure 14: Plot of ray curves for the plane oriented N65°.

Figure 15 :

 15 Figure15: Inversed velocity curves of P-waves as a function of the radial distance from the borehole wall for all of the measurements and in the direction of the minor and major horizontal stresses.

Figure 16 :

 16 Figure 16: Image of the inversed velocity field.

Figure 17 :

 17 Figure 17: Photograph taken from the head of the borehole.

Figure 1 :Figure 2 :

 12 Figure 1: Schematic diagram of the measurement method developed.

Figure 4 :

 4 Figure 4: P-wave velocity as a function of the radial distance for the numerical model (blue triangles) and the manufactured multilayer test bed (red crosses).

Figure 5 :

 5 Figure 5: Superimposed layers of the test bed and measurement device on the surface.

Figure 8 :

 8 Figure 8: Schematic diagram of the UPTO ultrasonic probe developed.

Figure 9 :Figure 11 :Figure 12 :Figure 13 :

 9111213 Figure 9: Meuse / Haute-Marne URL and experimental borehole location (courtesy of Andra (on left); Installation of the borehole probe (on right).

Figure 14 :Figure 15 :Figure 16 :

 141516 Figure 14: Plot of ray curves for the plane oriented N65°.

Figure 17 :

 17 Figure 17: Photograph taken from the head of the borehole.

Freewheeling probe head 14. Centralizers 15. Probe body 16. Power supply and transfer cables

  Figure 7: Comparison between the apparent velocities, expected and resulting from the inversion.

		5500 5500	16
	P-wave velocity (m/s) A. B. C.	3500 4000 4500 5000 3500 4000 4500 5000 1 Fixed part 14 Mobile Rotation part 2 A Mobile Translation part	3 4 5 1. Amplification Control 5 3 4 6 6 2. Magnetometer 3. Motors 4. Reduction gears 1. Amplification Control 2. Magnetometer 3. Motors 4. Reduction gears	7	9	10 12 9	MANUFACTURED (wanted) INVERTED MEASURED (apparent) 11 8 12 15 9 C B 9. Pneumatic claws 10. Ultrasonic transmitter 11. Ultrasonic receiver 12. Chocks 9. Pneumatic claws 10. Ultrasonic transmitter 11. Ultrasonic receiver 12. Chocks	13	9
			0.01 0.01	5. 5.	0.02 0.02 Encoders Encoders	0.03 0.03	0.04 0.04		0.05 0.05 13. Freewheeling probe head 13.
				6. 6.	Depth (m) Gears Gears				14. Centralizers
				7. 7.	Worm screw (left hand thread) Worm screw (left hand thread)	15. Probe body
				8. 8.	Worm screw (right hand thread ) Worm screw (right hand thread )	16. Power supply and transfer
										cables
	Figure 6: P-wave apparent velocity curve as a function of the inter-sensor distances for the
					multilayer test bed.

ACKNOWLEDGEMENTS

This work was supported by the European Union (EU) as part of the NF-PRO FI6W-CT-2003-02389 project (co-funded by the European Commission under the Euratom Research and Training Programme on Nuclear Energy within the Sixth Framework Programme: 2002-2006). The authors thank ANDRA (French National Radioactive Waste Management Agency) for its help with the experiment performed at the Meuse/Haute-Marne (Bure) Underground Research Laboratory.

TABLE OF FIGURES

Figure 10: Layout of measurements around the borehole (the radial lay out of the measurements is not realistic, it corresponds exclusively to the projected inter-sensor distances) with, in red, the discarded measurements.


Ghislain Gassier
email: ghislain.gassier@im2np.fr

Francoise Briolle
email: briolle@cpt.univ-mrs.fr

Gilles Chabriel
email: gilles.chabriel@im2np.fr

Jean Barrère
email: jean.barrere@im2np.fr

Passive Covert Radars using CP-OFDM SFN. Reference signal recovery from blind beamforming

A passive Coherent Location (PCL) system uses the signal transmitted by so called illuminators-of-opportunity in the environment for the illumination of targets. It is then necessary to recover the original transmitted signal to be compared to the targets echoes. With CP-OFDM transmissions, it is quite easy to recover the original data and then to reconstruct the original transmission. But in the case of a Single Frequency Network , characterized by the presence of several transmitters using the same carrier frequency to broadcast the same signal, it is necessary to use directional sensors or spatial filter. In this paper, we propose two blind solutions : classical beamforming and CAPON filtering. Results of real measurements are presented.

I. INTRODUCTION

A passive Coherent Location (PCL) system or Passive Covert Radar is an emerging and highly promising technology, which may be used in many application areas in the future [START_REF] Kuschel | Passive radar from history to future[END_REF], [START_REF] Thomas | HF passive bistatic radar potential and applications for remote sensing[END_REF], [START_REF] Song | Overview of passive radar technology based on non-cooperative radar illuminator[END_REF], [START_REF] Griffiths | Passive coherent location radar systems. Part 1: Performance prediction in Radar, Sonar and Navigation[END_REF].

A passive radar exploits the availability of so called illuminators-of-opportunity in the environment for the illumination of targets. The illuminators, normally used for radio or TV broadcasting, GSM transmissions, Wifi, etc. could be associated to a Passive Covert Radar [START_REF] Falcone | Experimental results for OFDM WiFi-based passive bistatic radar in Radar Conference[END_REF], [START_REF] Sun | Applications of passive surveillance radar system using cell phone base station illuminators[END_REF], [START_REF] Neyt | Feasibility of STAP for passive GSM-based radar in Radar[END_REF]. The keybenefits of passive radar are, apart from its covert operation, the low-cost nature due to saving on expensive transmitters as well as detection of low-altitude and/or stealthy targets.

A. Context

A typical Passive Covert bi-static Radar scene is illustrated in Fig. 1, where a standard transmitter illuminates a scene composed of different static reflectors due to topography, to buildings, . . . , and one or several moving targets (here an aircraft). A receiving antenna collects the radiations from the transmitter (direct path) and from the different reflectors.

Usually, a passive radar is equipped with two receiving channels: reference and surveillance channels. The signal from the reference channel, is used as the original transmitted signal (direct path) and compared to the surveillance signal.

Nowadays broadcast transmissions are changing more and more from analogue to digital. With CP-OFDM signals, used for Digital Audio and Video Broadcasting , Wifi transmissions, etc., it is possible to recover the transmitted digital data and then, to synthesize a sampled version of the original transmitted signal. The main difficulty of such a system will be first, to extract from the surveillance signal the only direct path and second, the very weak targets echoes. The direct path, used as a reference signal, will be compared to the target echoes in order to estimate the delay and the doppler effect.

After a brief description of the structure and properties of the CP-OFDM signals (SEC II), we will recall the main steps to process the demodulation of the signal in order to recover the original transmitted data (SEC III). The single-frequency network (SFN) broadcast is characterized by the presence of two or more transmitters using the same spectrum to broadcast the same signal. It is then necessary to use spatial filtering to recover the original data. Two spatial filters will be elaborated and results on real data will be presented in SEC IV.

II. CP-OFDM SIGNALS

The analog CP-OFDM signal is a sequence of OFDM blocks of duration T s ; the structure of each OFDM block relies on an orthogonal frequency division multiplexing (OFDM), which consists in a modulation method of orthogonal sub-carriers.

Each OFDM block is composed of two parts : the useful part of duration T u and a redundant part of duration T cp , called the cyclic prefix (CP).

The useful part of the OFDM block #m is composed of K orthogonal sub-carriers f k carrying a complex data d m k which can only take M different numerical complex values, depending on the M -array QAM modulation used. Note that orthogonality (on the duration T u ) between any pair of sub-carriers requires that f k = k Tu .

The cyclic prefix (CP) is obtained by copying the last part of the block and attaching it in front of the useful part. The duration T cp of the cyclic prefix is usually assumed to be greater than or equal to the delays of the echoes ; then, it allows to eliminate inter blocks (ISI) and inter carrier (ICI) interferences.

Assuming that the emitter time origin is taken at the emission of the first OFDM block (m = 0), we can express, in baseband, the m th emitted OFDM block as:

d m (t) = K-1 k=0 d m k exp(2iπ k Tu t) mT s < t ≤ (m + 1)T s .
And the direct path (the reference signal), sequence of blocks

d(t) = M m=0 d m (t),
Then, the signal d(t) is upconverted using the carrier frequency F 0 and the transmitter sends d(t)e 2iπF0t .

The main advantage of CP-OFDM communication technique is its ability to deal with multipath propagation without complex equalization filters. We will detail the process to recover the digital data {d m k } from the mixture received on the antenna of the passive radar.

III. CP-OFDM DEMODULATION

The signal from the surveillance channel is a mixture of the original transmitted signal (direct path) and echoes from static reflectors and moving targets, corrupted by sensor noise (Fig. 1). In baseband, this signal can be written as :

s surv (t) = A d (t -t 0 ) + i α i d (t -t 0 -τ i ) + j β j d (t -t 0 -τ j ) e i2πf d j t e 2iπF0t + η(t) (1)
where

• t 0 is the propagation delay between the receiving antenna and the transmiter,

• τ i + t 0 and τ j + t 0 are the delays between the receiving antenna and respectively the i th static reflector and the j th target,

• f d j is the frequency shift due to the radial velocity of the target j (Doppler effect),

• A, α i , β j are complex factors. α i and β j are related to positions and radar cross section (RCS) of the different reflectors.

• η(t) is the sensor noise contribution assumed to be centered.

Let denote d ref (t) the reference signal, proportional and delayed version of d(t), the surveillance signal can be written as:

s surv (t) = d ref (t) + i a i d ref (t -τ i ) + j b j d ref (t -τ j ) e i2πf d j t e 2iπF0t + η(t).

A. Process

The process to synthesize a sampled version of the reference signal d ref (t) performs the following 7 steps :

1) Baseband downconversion The surveillance signal s surv (t) is analogically multiplied by e -2iπF R 0 t to get the signal in the baseband:

x surv (t) = d ref (t) + i a i d ref (t -τ i ) + j b j d ref (t -τ j ) e i2πf d j t e 2iπf CF O t + η(t),
where f CF O := F R 0 -F 0 is the carrier frequency offset (CFO) between the transmitter and the receiver local oscillators. A method to estimate the CFO can be found e.g. in [START_REF] Sliskovic | Carrier and sampling frequency offset estimation and correction in multicarrier systems[END_REF].

2) Sampling

The duration of each OFDM block is T s = N s .∆, with ∆ the sampling period, given by the CP-OFDM standard. We consider here that the sampling frequency offset (SFO) effect is negligible. However, some estimation and correction methods can be found in [START_REF] Sliskovic | Carrier and sampling frequency offset estimation and correction in multicarrier systems[END_REF], [START_REF] Shim | OFDM carrier frequency offset estimation methods with improved performance[END_REF], [START_REF] You | An improved sampling frequency offset estimator for OFDM-based digital radio mondiale systems[END_REF]. If we assume that τ k = n k ∆ is an integer multiple of the sampling period, then the sampled surveillance signal yields

x surv (n) = d ref (n) + i a i d ref (n -n i ) + j b j d ref (n -n j ) e i2πf d j n e 2iπf CF O n + η(n).

3) Synchronization

The goal of this synchronization is to detect the beginning of a CP-OFDM block. Each CP-OFDM block is composed of K orthogonal sub-carriers. P of them, called pilots frequencies, carry a specific value {d p } . The values and the frequencies {f p } of the pilots are given by the DVB-T standard.

The synthesized sampled time sequence of pilots, over N OFDM blocks, is

x P (n) = P p=1
d p e 2iπfpn n = 1 . . . N.N s , Nota : The DVB-T norm imposes that the frequencies of some pilots (continual and TPS) are the same for all OFDM blocks. The frequency of scattered pilots differs from one block to another, with a periodicity of four blocks. Here, for simplicity, this detail is omitted (see [?]).

The position m 0 of the first peak of the intercorrelation sequence between the sampled surveillance signal x surv (n) and the sequence of pilots x P (n) gives the position of the first OFDM block.

Then, the surveillance signal can be written as a sequence of N OFDM blocks, with the block #m 4) Equalization

For each OFDM block, the cyclic prefix will be removed. The block #m is equal to

x m surv (n) = K-
d m k = dm k × H( k T u ). (2)
This equalization will be done for each block. Figures 6 or 7 show the real part versus the imaginary part (constellation diagram) of the equalized data {d k } for N OFDM blocks. We can note here that the data {d k } take only 16 different complex values plus 4 real values used for the pilots.

5) CFO correction

The continual and TPS pilots have the same frequency f p for each OFDM block. The Carrier Frequency Offset will be estimated by the phase shift of data pilots through N OFDM blocks. The average phase shift for 60 OFDM blocks is plotted on figure 3. A linear regression estimation will give the slope S of the line and the estimation of f CF O = S 2π . Multiplying the x surv (n) by e -2iπf CF O n will correct the CFO effect.

Then after correction, the equalization process will be done again as it is described above. 6) Mapping Fig. 3: Average of the phase shift fover 60 OFDM blocks

The binary flow was encoding with a 16-QAM. We have then to map the equalized set {d k } to 16 complex values plus 4 real ones to obtain the original transmitted data.

7) Reconstruction of the reference signal

From the original transmitted data set,the sampled useful part, of duration T u , is obtained by FFT for each OFDM block. Then by copying the last part and attaching it in front, we will reconstruct the OFDM block of duration T s .

IV. DEMODULATION IN SFN NETWORK CONFIGURATION: RESULTS ON REAL MEASURMENTS

A. Context

A single-frequency network (SFN) is a broadcast network characterized by the presence of several transmitters using the same carrier frequency to broadcast simultaneously a same signal.

The campaign of measurements, we present, uses three DVB-T SFN illuminators (8K-mode, 514 MHz). The receiver is an array composed of 16 sensors. The positions of the different illuminators and the receiver are described on figure 4. The direct path of the illuminator #1 is followed by the strong direct paths of the illuminators #2 then #3, which can be seen as "strong echoes".

B. Modelling

With P transmitters in the SFN, the surveillance signal, can be modelled as

s surv (t) = P p=1 A p d(t -t p ) + ( P p=1 i α p,i d(t -τ p,i ) + P p=1 j β p,j d(t -τ p,j t)e i2πf d p,j t e 2iπF0t + η(t)
where d(t -t p ) is the direct-path signal of p-th transmitter and A p the corresponding complex amplitude; α p,i and β p,j are the complex amplitude of static reflectors and the targets echoes and τ p,i and τ p,j , the delays.

Here, the signal of interest (the mos powerful) is the direct path from the nearest transmitter, called the reference signal. Then the surveillance signal, after baseband downconversion, can also be written as,

x surv (t) = d ref (t) + i,j a i d ref (t -τ i ) + d ref (t -τ j ) e i2πf d j t P p=2 A p d(t -t p ) + i,j α p,i d(t -τ p,i ) + β p,j d(t -τ p,j t)e i2πf d p,j t e 2iπf CF O t + η(t).
The main issue, when SFN illuminators are considered is that the P "strong echoes" due to these different transmitters produce a deep modulated channel frequency response so that the interpolation process fails.

Thanks to the CP-OFDM structure of the DVB-T signals, the synchronization can still be achieved by correlation with a pilot sequence. But if any spatial process is done to enhance the direct path from the nearest transmitter, then it will be impossible to recover the original data interpolating the pilotfree frequencies. The results (Fig. 5) are obtained when a single sensor is considered.

C. Spatial blind enhancement

The goal is here to enhance the signal emitted from the nearest transmitter.

Let us consider

• {x 1 surv (n), . . . x M surv (n)}, the surveillance signals recorded on the M sensors of the antenna,

• x s , the observations vector,

x s = [x 1 surv (n), . . . x M surv (n)] T , • R x spatial covariance matrix of the observations, R x = x s .x s H .
To enhance the signalfrom the direction of arrival (DOA) θ 1 corresponding to the nearest transmitter, we have to estimate where u(θ 1 ) = l sin(θ 1 )/λ, l is the distance between two adjacent sensors and λ the wavelength corresponding to the carrier frequency. We show in the full paper that this steering vector can be blindly (without any geometrical or DOA consideration) estimated using both the sequence of pilots x P (n) used for the synchronization and the synchronized surveillance signals x m surv (n). We show in particular that d(θ 1 ) can be estimated using

d(θ 1 ) = 1, Γ 2 (0) Γ 1 (0) . . . , Γ M (0) Γ 1 (0) T where Γ m (0) := E{x * P (n)x m surv (n)}.
Now, from this estimated steering vector corresponding to the nearest transmitter, one can implement a spatial filtering in order to produce a beam where the contributionof the strong echoes due the presence of the other transmitters is reduced. In this study, we consider two different spatial filters w(θ 1 ):

• The classical beamformer:

w CB (θ 1 ) = d(θ 1 )

• The CAPON beamformer:

w CAPON (θ 1 ) = R x -1 d(θ 1 ) d(θ 1 ) H R x -1 d(θ 1 )
The two following figures shows the constellation diagrams obtained with the different spatial filters before the mapping step. One can see that in particular with the CAPON filtering, the quality of the obtained constellation is sufficiently good so that the last mapping process conducts to a perfect reference signal recovery. These figures have to be compared with Figure 5 where no spatial processing was implemented. The Single Frequency Network broadcast is used in many countries : Germany, North of Europe, China, etc.. . For such a network, the reference signal recovery needs a spatial processing. Two solutions were implemented : classical beamforming and spatial CAPON filter with a blind estimation of the steering vector of interest. The results obtained with real data indicate that the CAPON filter gives better results; the constellation diagram of the estimated data exhibits a satisfactory variance. Now, the computation of such a filter requires many operations while the implementation of classical beamforming is quite simple.

The main advantages for the Passive Coherent Location using CP-OFDM in a SFN mode is to have the targets illuminated by several transmitters. It is then possible to localize accurately the targets.

Fig. 1 :

 1 Fig. 1: Passive radar scene.

Fig. 2 :

 2 Fig. 2: Structure of a CP-OFDM block.

p

  where dm k = H m ( k Tu )d m k is the complex amplitude of the k th sub-carrier. This value can be easily obtained by Inverse Fast Fourier Transform. The values {d m p } of the pilots are known, and the frequency response of the propagation channel is known for the pilots frequencies {f m p }The frequency response of the propagation channel, H( k Tu ) for k = p, is obtained by linear interpolation. Then, an estimation of the transmitted data can be obtained by equalization,

Fig. 4 :

 4 Fig. 4: Positions of the illuminators and the PCL receiver.

Fig. 5 :

 5 Fig. 5: Constellation diagram of the equalized data (without any spatial processing)

Fig. 6 :Fig. 7 :

 67 Fig. 6: CB-equalized data constellation diagram (0.4 Hz CFO compensation)

ACKNOWLEDGMENT

The campaign of measurements was achieved by Thales Air Systems, ONERA, DGA (Direction Générale de l'Armement) and CEAM (Centre d'Expériences Aériennes Militaires) jointly with the BWB (Bundesamt fr Wehrtechnik und Beschaffung) and Fraunhofer Institute. We are particularly grateful to J. Colin, from the DGA, for giving us access to the data.


Gilles Chabriel
email: gilles.chabriel@im2np.fr

Jean Barrère
email: jean.barrere@im2np.fr

Ghislain Gassier
email: ghislain.gassier@im2np.fr

Francoise Briolle
email: briolle@cpt.univ-mrs.fr

Passive Covert Radars using CP-OFDM signals. A new efficient method to extract targets echoes

published or not. The documents may come

Passive Covert Radars using CP-OFDM signals. A new efficient method to extract targets echoes

I. INTRODUCTION

A Passive Coherent Location (PCL) system or Passive Covert Radar is an emerging and highly promising technology, which may be used in many application areas in the future [START_REF] Kuschel | Passive radar from history to future[END_REF], [START_REF] Thomas | HF passive bistatic radar potential and applications for remote sensing[END_REF], [START_REF] Song | Overview of passive radar technology based on non-cooperative radar illuminator[END_REF], [START_REF] Griffiths | Passive coherent location radar systems. Part 1: Performance prediction in Radar, Sonar and Navigation[END_REF].

A passive radar exploits one (bistatic case) or several (multistatic case) so called "illuminators-of-opportunity" in the environment to detect moving targets through their Doppler . The illuminators, normally used for radio or TV broadcasting, GSM transmissions, Wifi, etc. could be associated to a Passive Covert Radar [START_REF] Falcone | Experimental results for OFDM WiFi-based passive bistatic radar in Radar Conference[END_REF], [START_REF] Sun | Applications of passive surveillance radar system using cell phone base station illuminators[END_REF], [START_REF] Neyt | Feasibility of STAP for passive GSM-based radar in Radar[END_REF].

The main advantages of such systems, with respect to conventional (active) radars, include low cost architecture, low energy requirement, detection of stealth targets and potentially null probability of interception.

Besides these incontestable practical advantages, a major problem is that Doppler echoes are in general more difficult to detect than with an active (monostatic) radar, since the bistatic angle is in general greater than a mono-static angle, and the target radar cross section (RCS) lower. Continual presence of the line of sight (LOS) emission with several orders of magnitude greater than the Doppler echoes may also bury the echoes of interest.

Reduction of the unwanted zero-Doppler echoes (defined here as "clutter") is then particularly important in such systems if one wants to expect an efficient target detection.

Nowadays broadcast transmissions are changing more and more from analog to digital. With CP-OFDM signals, widely used for communications transmissions as Digital Audio and Video Broadcasting, Wifi transmissions, etc., it is possible to recover the original transmitted digital data and then, to synthesize a sampled version of the original transmitted signal. This signal will be used as a reference signal.

The main difficulty of PCL is to extract from the surveillance signal the very weak target's echoes which will be compared to the reference signal (direct path) in order to estimate the delay and the doppler effect.

After a presentation of Passive Covert Radar using CP-OFDM and a brief presentation of the process to recover the original transmitted signal (SEC:II), a new efficient method to extract targets echoes from the clutter (SEC:III) and results are presented (SEC:IV) on simulated data and on real data.

II. POSITION OF THE PROBLEM

A. Modelling

A typical Passive Covert bi-static Radar scene is illustrated in Fig. 1, where a standard transmitter illuminates a scene composed of different static reflectors due to topography, buildings, . . . , and one or several moving reflectors of interest or targets (here an aircraft). A receiving antenna collects the radiations from the CP-OFDM transmitter and from the different reflectors. Usually, a passive radar is equipped with two receiving channels: reference and surveillance channels. The signal from the reference channel, is used as the original transmitted signal (direct path) and compared to the surveillance signal.

On the receiving antenna, we consider the baseband received signal s(t) as a linear combination of the CP-OFDM emitted signal s e (t) following the direct path, and of its static and Doppler echoes corrupted by sensor noise.

We can write: [START_REF] Kuschel | Passive radar from history to future[END_REF] where:

s (t) =αs e (t -t 0 ) + i β i s e (t -τ i -t 0 ) + j γ j s e (t -τ j -t 0 ) exp i2πf dj (t) + η(t),
• t 0 is the propagation delay between the receiving antenna and the transmitter,

• τ i + t 0 and τ j + t 0 are the delays between the receiving antenna and respectively the i th static reflector and the j th target,

• f dj is the frequency shift due to the radial velocity of the target j (Doppler effect),

• α, β i , γ j are complex factors. β i and γ j are related to positions and RCS of the different reflectors.

• η(t) is the sensor noise contribution assumed to be centered.

NOTA: the model of the baseband signal do not take into account the carrier frequency offset (CFO) effect due to the difference between the transmitter and receiver local oscillators. Different methods can be used to correct this effect ([10], [START_REF] Shim | OFDM carrier frequency offset estimation methods with improved performance, Broadcasting[END_REF], [START_REF] You | An improved sampling frequency offset estimator for OFDM-based digital radio mondiale systems[END_REF]). Unlike conventional active radar systems, the original emitted signal s e (t) is not directly available. Nevertheless, with a CP-OFDM transmission, one can recover perfectly the transmitted data and consequently synthesize the baseband version of the emitted signal. This signal will be called the reference signal denoted by s ref (t). This process is detailed in II-B.

Then, the received signal can be written as

s (t) =αs ref (t) + i β i s ref (t -τ i ) + j γ j s ref (t -τ j ) exp i2πf dj t + η(t).
In practice, the target's radar cross section (RCS) is very weak and consequently |γ j | is very small compared to |α| and of most of

|β i |.
It is then necessary to enhance, from the received signal s(t), the target's echoes signal

s tar (t) = j γ j s ref (t -τ j ) exp(i2πf dj t)
or equivalently to reduce the unwanted echoes, denoted clutter

s clut (t) = αs ref (t) + i β i s ref (t -τ i )
by a proper signal processing.

Finally, this ideally "clutter-cleaned" signal will be "compared" to the reference signal s ref (t) (ambiguity function) in order to estimate the radial speed and the bi-static distance of the different targets. Different methods have been proposed [START_REF] Colone | A multistage processing algorithm for disturbance removal and target detection in passive bistatic radar Aerospace and Electronic Systems[END_REF], [START_REF] Poullin | Passive Radar using COFDM (DAB or DVB-T) Broadcasters using as Opportunistic Illuminators[END_REF]. The new efficient proposed method here will be based on the property of CP-OFDM communications that is a precise estimation of the propagation channel which allow echoes cancellation.

B. Construction of the reference signal

The structure of CP-OFDM signal allows a simplified correction of the multipath effect in order to recover the original transmitted data. Indeed, with such a signal, it is possible to get an accurate estimation of the propagation channel and then, by inverse filtering, to get the original data. These data will be used to synthesize a sampled version of the original signal, used as a reference signal. Let us briefly recall the relevant characteristics of a CP-OFDM emission.

1) CP-OFDM signals:

CP-OFDM signal relies on an orthogonal frequency division multiplexing (OFDM), which consists of a modulation method of closely spaced orthogonal sub-carriers. The OFDM analog signal is a sequence of symbols. Each symbol, of duration T s , is composed of K orthogonal sub-carriers ν k carrying the complex data d m k . A data can only take M different numerical complex values depending on the M -array QAM modulation used. A symbol consists of two parts: a useful part of duration T u (carrying the data) and a redundant part of duration ∆ so-called the cyclic prefix (CP). The cyclic prefix (CP) is obtained by copying the last part of the symbol and attaching it in front of the useful part. CP allows to eliminate inter symbol (ISI) and inter carrier (ICI) interferences induced by significant echoes of delay less than ∆. Note that orthogonality (on the duration T u ) between any pair of sub-carriers requires that ν k = k Tu .

Assuming that the transmitter time origin is taken at the emission of the first symbol (m = 0), we can express the useful part of the m th emitted symbol as:

s m e (t) = K-1 k=0 d m k exp(2iπ k Tu t) mT s + ∆ < t ≤ (m + 1)T s .
2) Estimation of the propagation channel: to simplify the exposure, let us initially consider a noise-free case without echoes from mobile targets. The received signal is then the sum of the direct path and echoes from static reflectors. After synchronization onto the first emitted symbol and CP elimination, the received signal corresponding the m th emitted symbol can be written as a function of time, for the duration T u , as

s m (t) = K-1 k=0 dm k exp(2iπ k Tu t), mT s + ∆ < t ≤ (m + 1)T s ,
where dm k = H m (ν k )d m k is the transmitted data (complex amplitude of the k th sub-carrier). The unknown coefficient H m (ν k ) is the frequency response of the propagation channel at the frequency ν k at the time mT s . Note that under a block-fading channel assumption, one can consider that H m (ν k ) remains constant over several symbols.

Let us sample the signal s m (t) at the frequency K Tu , given by the CP-OFDM standard; we obtain the following numerical sequence:

s m [n] = K-1 k=0 dm k exp(2iπ kn K ), 0 ≤ n ≤ K -1.
The coefficients dm k are then obtained by discrete Fourier transform through a fast Fourier transform algorithm (FFT):

dm k = N -1 n=0 s m [n] exp(-2iπ kn K ), 0 ≤ k ≤ K -1.
Most of CP-OFDM standard, like DVB, requires that, inside a symbol m, P particular sub-carriers ν p called pilot frequencies, carry the fixed data {d m p } p see e.g. [13]. Now, assuming that the coefficients {d m p } p,m are known, one can compute the frequency response {H m (ν p )} p,m of the propagation channel at the frequencies {ν p } p,m :

H m (ν p ) = dm p d m p , p ∈ 1, • • • , P.
The frequency response of the channel for the others frequencies will be estimated by interpolation and the estimation of the original data can be processed by inverse filtering,

d m k = dm k H m (ν k ) , k ∈ 0, • • • , K -1.
Actually, both the presence of noise η(t) and Doppler echoes s tar (t) introduce errors during this estimation process. However, when the signal-to-noise ratio is sufficient and when the levels |γ i | of the Doppler echoes are weak enough, mapping the estimated data in accordance with the QAM constellation allows to correct most of these errors. We will consider here that after this process, the data are perfectly recovered.

We can already note that if the original transmitted data are perfectly recovered, then the propagation channel can be estimated with a very strong accuracy, for each symbol.

3) The reference signal: applying an inverse fast Fourier Transform (IFFT) onto the sequence {d

s ref = {s 0 ref , s 1 ref , • • • , s m ref , • • • }.

III. CLUTTER REMOVAL METHOD

The method will rely on the hypothesis that the effect due to the static reflectors (including the direct path), will be the "same over several symbols (the duration T s of a symbol is about 1ms).

Since, with CP-OFDM signals, it is possible to recover the original data with a great accuracy, it is also possible to estimate, with a great accuracy, the frequency response H m (ν k ) of the propagation channel for each symbol m. If it is possible to estimate the propagation channel H c (ν k ) of the clutter ( echoes from static reflectors and direct path), then it will be easy to synthesize the clutter:

S clut (ν k ) = H c (ν k ).S ref (ν k ).
Then, the clutter will be removed from the surveillance signal to obtain only the target's echoes.

A. Method

from the general noisy case (II-A) involving clutter and Doppler echoes, the sampled version of the baseband received signal s(t) is:

{s} = {s clut } + {s tar } + {η}.
Now, considering any symbol m without its CP, we have

s m [n] = s m clut [n] + s m tar [n] + η m [n], 0 ≤ n ≤ K -1.
Under the block-fading channel assumption, s m clut [n] is expressed as

s m clut [n] = K-1 k=0 H c (ν k )d m k exp(2iπν k )), 0 ≤ n ≤ K -1,
where H c (ν k ) is the channel frequency response due to the clutter, assumed to be constant for the M symbols.

Computing the FFT of the m th received symbol

s m [n], one obtains S m [k] = dm k = H m (ν k )d m k
where dm k are the received data and d m k are the original one. H m (ν k ) is the frequency response of the propagation channel for the symbol m, given by

H m (ν k ) = dm k d m k
The FFT of the m-th received symbol can also be written, for We can divide the previous equation (Eq.2) by d m k to get

0 ≤ k ≤ K -1, as S m [k] = H m (ν k )d m k = H c ( k Tu )d m k + S m tar [k] + N m [k], (2)
Ĥm (ν k ) = H c (ν k ) + S m tar [k] d m k + N m [k] d m k , 0 ≤ k ≤ K -1
One can demonstrate that

Ĥc (ν k ) = 1 M M -1 m=0 Ĥm (ν k )
is an unbiased estimator of the channel frequency response H c (ν k ) due to the clutter effect.

From Ĥc (ν k ) one can obtain an estimation of S m tar [k] by computing

Ŝm tar [k] = S m [k] -Ĥc (ν k ).S ref [k].

B. Process

The different stages of the proposed frequency clutter removal method can be found in Fig. 2. IV. RESULTS

A. Simulated data

We use 90 symbols of an original CP-OFDM signal, encoded with the DVB-T standard (8K, QAM-64, T ∆ = T u /16). The sensors noise is a white gaussian noise and the signal to noise ratio is equal to 40 dB.

The echoes from 100 static reflectors are added to the original CP-OFDM signal, delayed randomly by τ i and amplified by β i such as τ i < T ∆ and 20log 10 ( α βi ) < -30dB, with α the amplitude of the original signal. The echoes direction (phase) is random.

The echoes of 3 targets are added to this mixture: Target #1 : distance to the receptor 10 km, doppler effect 20 Hz and α A1 =-50 dB, with A 1 the echo's amplitude, Target #2 : distance to the receptor 12,5 km, doppler effect 125 Hz and α A2 ==-50 dB, with A 2 the echo's amplitude, Target #3 : distance to the receptor 30 km, doppler effect -150 Hz and α A3 ==-60 dB, with A 3 the echo's amplitude. Two of the targets are very close to the PCL receiver (12 km and 12.5 km). One of them moves very slowly, with a doppler echo equal to 20 Hz. The third target is far from the receiver (30 km) and the attenuation of the echo is 60 dB.

The

1) Results: the representation of the ambiguity function, [START_REF] Woodward | Probability and information theory, with applications to radar[END_REF], [START_REF] Chevalier | Principles of radar and sonar signal processing[END_REF], defined as

χ(x, ν) = +∞ -∞ s ref (t)s * tar (t -x/c)e 2iπνt dt
will give an image ( doppler frequency versus bi-static distance) where the maxima of the ambiguity function the target will be plotted in red. Two methods were used : our new method (presented in SEC:III) and the least squares based method proposed by Colone et al. [START_REF] Colone | A multistage processing algorithm for disturbance removal and target detection in passive bistatic radar Aerospace and Electronic Systems[END_REF], called ECA-batches, where the signal was cut in thirty parts of three symbols each. The results are presented in fig. 3 and4.

The results obtained with the new efficient method allow to detect the three targets on the simulated data, even when the target is moving slowly (target #2). Withthe ECA-batches method it was impossible to detect this target, which was seen as a clutter (zero doppler) echo.

B. Real data

The campaign of measurements was achieved using illuminators in a Single Frequency Network. The broadcast signal, encoded with DVB-T standard (8K, QAM-16, T ∆ = T u /16, is modulated at the frequency of 514 MHz. The receiver was an array composed of 16 sensors located close to an airport and an highway.

The cross-ambiguity function (extended to the multi-sensor case) presented in figure 5 is obtained from 82 OFDM symbols. The method presented here seems to be promising. First, because its implementation is very simple and requires few computation. The results on simulated data show that the method is reliableand results on real data are hopeful.

Fig. 1 :

 1 Fig. 1: Passive radar scene.

  m k } k to obtain the useful part of the symbol and attaching in front the CP insertion, one obtains the sampled version {s m ref [n]} n of the original transmitted CP-OFDM symbol m. A sequence of symbols {s m ref [n]} n,m will provides the whole reference signal {s ref } we seek:

  where S m tar [k] and η m [k] are respectively the discrete Fourier transform of s m tar [n] and η m [n].

Fig. 2 :

 2 Fig. 2: Frequency clutter removal method.

Fig. 3 :

 3 Fig. 3: New efficient method -ambiguity function

Fig. 5 :Fig. 6 :

 56 Fig. 5: New efficient method -ambiguity function

  parameters of the simulation are presented in the table below.

		Dop. freq. (Hz) At. (dB) Range (km)
	Echoes (100)	0	>30	< 16
	Target #1	20	50	10
	Target #2	125	50	12,5
	Target #3	-150	60

ACKNOWLEDGMENT

The campaign of measurements was achieved by Thales Air Systems, ONERA, DGA (Direction générale de l'armement) and CEAM (Centre d'Expériences Aériennes Militaires) jointly with the BWB (Bundesamt fur Wehrtechnik und Beschaffung) and Fraunhofer Institute. We are particularly grateful to J. Colin, from the DGA, for giving us access to the data.