Jean Auriol

Roman Shor

Silviu-Iulian Niculescu

Nasser Kazemi

Estimating Drill String Friction Parameters: Comparing Performance of Model-Based Estimators to a Data-Driven Neural Network

In this paper, we consider the torsional motion of a drilling system and propose three algorithms to estimate the friction factors that characterize the interaction between the drill pipe and the wellbore walls (Coulomb source terms). This is essential to design the next generation of stick-slip mitigation controllers, to develop real-time wellbore monitoring tools, and to enable effective toolface control for directional drilling. We propose two model-based algorithms (an adaptive observer and a recursive dynamics framework) and a machine learning-based algorithm to estimate friction parameters, all of them presenting advantages and drawbacks. The performances of our modelbased estimators are finally compared with this data-driven neural network.

I. INTRODUCTION

Extraction of resources in the earth's subsurface -oil, gas, minerals, and thermal energy -necessitates the drilling of long slender boreholes from the surface to the subsurface target. The diameters of these wells range from 10 to 50 cm, and lengths can frequently exceed 10,000 m, leading to mechanical systems with extreme aspect ratios. These drill strings consist of sections of steel drill pipe and stiffer drill collars. Transfer of energy (rotational and axial) often leads to complex dynamic behaviors, which often manifest themselves as drill string vibrations and has been extensively studied [START_REF] Aarsnes | Axial and torsional self-excited vibrations of a distributed drill-string[END_REF], [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF] or [START_REF] Saldivar | A control oriented guided tour in oilwell drilling vibration modeling[END_REF] for a more complete review.

A particular vibration pattern that is prevalent is called stick-slip, where self-excited torsional oscillations manifest themselves in the drill string and are caused by friction along the wellbore between the drill pipe and the wellbore wall and by the bit-rock interaction. These stick-slip oscillations are generally characterized by a series of stopping -"sticking" -and releasing -"slipping" -events of the bit.

This phenomenon has been studied extensively and has been explained through a non-linear, or velocity-weakening -effect of bit-rock interaction [START_REF] Brett | Uses and Limitations of Drillstring Tension and Torque Models for Monitoring Hole Conditions[END_REF], [START_REF] Halsey | Drillstring Vibrations: Comparison Between Theory and Experiments on a Full-Scale Research Drilling Rig[END_REF], through distributed friction between the drill pipe and the wellbore wall [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] or through a combination of both [START_REF] Auriol | Self-tuning torsional drilling model for real-time applications[END_REF], [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF], [START_REF] Richard | Self-excited stick-slip oscillations of drill bits[END_REF]. This is of particular importance in modern wellbores, which are rarely straight and must follow pre-planned well plans, ranging from simpler horizontal or deviated wells to complex threedimensional paths, thus increasing the effect of torque and drag.

During directional drilling operations, downhole measurements are usually available. They are obtained by combining accelerometers or gyroscopes and one to three-axis magnetometers. These sensors are sampled at frequencies between 1-100 Hertz by the downhole tool but only averaged or windowed values are transmitted to surface [START_REF] De Almeida Jr | A review of telemetry data transmission in unconventional petroleum environments focused on information density and reliability[END_REF]. Thus, these measurements are transmitted to the surface 1 Jean Auriol, and Silviu-Iulian Niculescu are with Université Paris-Saclay, CNRS, CentraleSupélec, Inria, Laboratoire des Signaux et Systèmes, 91190, Gif-sur-Yvette, France. 2 Roman Shor and Nasser Kazemi are with the University of Calgary, Calgary, Alberta, Canada. Corresponding author: jean.auriol@centralesupelec.fr infrequently (once per minute to once per hour) and with a significant delay (on the order of seconds). Although this is not a problem for human-in-the-loop operations, it becomes insufficient for automated solutions. Such delays and sampled measurements imply significant performance degradation when using feedback controllers that aim to compensate stick-slip oscillations [START_REF] Aarsnes | Benchmarking of industrial stick-slip mitigation controllers[END_REF]. Real-time estimation of friction parameters [START_REF] De Wit | A new model for control of systems with friction[END_REF] enable feedforward [START_REF] Aarsnes | Avoiding stick slip vibrations in drilling through startup trajectory design[END_REF] control as well as a new generation of controllers based on the backstepping methodology [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] and thereby allow better performance.

In the present work, we propose three strategies (two of them being completely original) to estimate these Coulomb friction factors. Two of them are model-based (adaptive observer and recursive dynamics framework), while the last one is based on machine-learning. The model we use for the two first approaches is adjusted from [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF]: the bit is off-bottom (we do not consider bit-rock interaction), and Coulomb friction is included as a source term implemented as an inclusion.

This paper is structured as follows: it begins with a brief derivation of the drill string model in section II, followed by the description of three methods to estimate friction: 1) Section III describes the adaptive observer, 2) Section IV describes the recursive dynamics framework, and 3) Section V describes the machine learning approach Finally, in Section VI the three approaches are compared with data generated by the model described in section II.

II. MODEL UNDER CONSIDERATION

In this section, we present the mathematical model developed in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] to describe the torsional dynamics of the drill string. The proposed model is relatively simple (which makes it usable for control and estimation perspectives) and realistic since it has been validated against field data. The main assumptions we use are the following:

• The torsional motion of the drill string is the dominating dynamic behavior. • Uniform axial motion. No distributed axial dynamics.

• We consider that the bit is off-bottom, i.e., there is no bit-rock interaction. • The transition from static to dynamic Coulomb friction is modeled as a jump, i.e., the Stribeck curve is assumed negligible.

• The effects of along-string cuttings distribution on the friction is assumed to be homogeneous. • The effect of the pressure differential, inside and outside the drill string, on the bending moment, is considered to be negligible.

ω TD τ(t,x) ω(t,x) x INC x=L S(ω) ω c ω r o (x)μ s F N (x) 1 k t r o (x)μ k F N (x)

A. Torsional dynamics of the drill string

We consider that the torsional motion of the drilling device corresponds to the dominating dynamic behavior. The distributed model we present here is adjusted from [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] (where the reader can find the complete model derivation). This kind of representation is popular in the literature [START_REF] Aarsnes | Axial and torsional self-excited vibrations of a distributed drill-string[END_REF], [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF]. We give a schematic representation of the drill string in Fig. 1. Let us denote the angular velocity and torque as ω(t, x), τ (t, x), respectively, with (t, x) ∈ [0, ∞) × [0, L] (L being the length of the drill string). We have

∂τ (t, x) ∂t + JG ∂ω(t, x) ∂x = 0 (1) Jρ ∂ω(t, x) ∂t + ∂τ (t, x) ∂x = S(t, x), (2)
where, J is the polar moment of inertia, G the shear modulus and ρ the drill string density. The source term S is given as

S(t, x) = -k t ρJω(t, x) -F(t, x), (3)
where the damping constant k t is the viscous shear stresses and where F is a differential inclusion that corresponds to the Coulomb friction between the drill string and the borehole, also known as the side force. We use the following inclusion to implement the side force

   F(t, x) = r o (x)µ k F N (x), ω(t, x) > ω c , F(t, x) ∈ ±r o (x)µ s F N (x), |ω(t, x)| < ω c , F(t, x) = -r o (x)µ k F N (x), ω(t, x) < -ω c , (4)
where µ s is the static friction coefficient (i.e. the friction between two or more solid objects that are not moving relative to each other) and µ k is the kinetic friction coefficient (also known as dynamic friction or sliding friction, which occurs when two objects are moving relative to each other and rub together), ω c is the threshold on the torsional velocity where the Coulomb friction transits from static to dynamic. The function F N is the normal force acting between the drill string and the borehole wall. The function F(ω) ∈ ±r o (x)µ s F N (x) denotes the inclusion where

F(t, x) = - ∂τ (t, x) ∂x -k t ρJω(t, x) ∈ [-r o (x)µ s F N (x), r o (x)µ s F N (x)], (5)
and take the boundary values ±µ s F N (x) if this relation does not hold. The shape of the friction source term is illustrated in Fig. 1.

The method to calculate normal force, to account for discontinuities of a multiple sectioned drill string and the Riemann invariants may be found in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF].

B. Boundary condition

We consider that the drilling device is actuated at the topdrive by a motor torque: τ m . To reach a given reference velocity set-point ω SP , we consider a standard PI control law [START_REF] Åström | Feedback systems: an introduction for scientists and engineers[END_REF]:

τ m = k p (ω SP -ω 0 ) + k i t 0 (ω SP (ξ) -ω 0 (ξ))dξ, (6)
where k p is a proportional gain and k i an integral gain. We denote J T D the topdrive inertia and where ω 0 verifies the following equation

∂ω 0 ∂t = 1 J T D τ m -τ (t, 0) . (7)
Finally, since the bit is off-bottom, we have τ (t, L) = 0.

C. Estimation problem

As detailed in the introduction, the design of the next generation of stick-slip mitigation controllers requires reliable estimations of the friction parameters. In the rest of the paper, we present three methods to estimate all (or part) of the friction parameters. We emphasize their advantages and their drawbacks. Data measurements correspond to the topdrive torque and velocity.

III. FIRST APPROACH: ADAPTIVE OBSERVER

The first approach we propose is adjusted from the adaptive state-observer proposed in [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF]. Although the results given in [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF] correspond to the case of a two-sections drilling device (i.e., the drill string is made of only two sections: the pipe and the collar), they can be extended to the case of a multisectional drill string, either by adjusting the theoretical approach given [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] or the new recursive framework proposed in [START_REF] Redaud | Output-feedback control of an underactuated network of interconnected hyperbolic pde-ode systems[END_REF]. This observer combines measurements from physical sensors (top-drive angular velocity ω 0 ) with the proposed model of the system dynamics. It provides reliable estimates of the torque and RPM states and of the side forces friction parameters when the bit is off bottom [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF]. We recall here the main ideas of this observer.

Let us denote with the • superscript the estimated states and e = ω0 -ω 0 the measured estimation error of the topdrive angular velocity. The observer equations given in [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF] in terms of Riemann invariants read as follows (the index i corresponds to the considered section of the drill-string) [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF].

ω0 = a 0 βp (t, 0) -ω0 + 1 I T D τ m -p 0 e, (8)
∂ αi ∂t (t, x) + (c t ) i ∂ αi ∂x (t, x) = Ŝi (t, x) -p i α (x)e, (9)

Requires the knowledge of ωc

Convergence guaranteed in the absence of friction terms [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] No proof of convergence for the adaptive part Real-time estimation robust to uncertainties and delays Persistent excitation for the adaptive part

∂ βi ∂t (t, x) -(c t ) i ∂ βi ∂x (t, x) = Ŝi (t, x) -p i β (x)e, (10)
The source term in each section is computed from the estimated states and friction factor:

Ŝi (t, x) = k t (α i (t, x) + βi (t, x)) + 1 J i ρ i Fi (t, x), (11)
where F has an expression identical to (4), the different variables being replaced by their estimates. The different boundary conditions now express as

α(t, 0) = 2ω 0 (t) -β(t, 0) -P 0 e, (12)
βi (t, x i+1 ) = αi (t, x i+1 )(1-Zi ) + 2 Zi βi+1 (t, x i+1 ) 1 + Zi -P 1 e, (13)
αi+1 (t, x i+1 ) = 2 Zi αi (t, x i+1 ) -βi+1 (t, x i+1 )(1-Zi ) 1 + Zi , (14)
β(t, L) = α(t, L), (15)
and the estimates of the friction factor are updated according to

μs (t) = -l s e, |min x ω(t, x)| ≤ ω c , 0, |min x ω(t, x)| > ω c , (16)
μk (t) = 0, |min x ω(t, x)| ≤ ω c , l k e, |min x ω(t, x)| > ω c , (17)
Finally, we use a saturation to improve the robustness of the approach: μs = max(μ s , μk ). The different constants and observer gains a 0 , p i α , p i β , p 0 , p 1 , P 0 , P 1 , l s , l k can be found in [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF]. The estimation procedure is summarized by Algorithm 1 Algorithm 1 Estimation of µ k and µ k with an adaptive observer. Require: τ m , ω 0 , ω c 1: Simulate system (8)-( 15) with the adaptive laws ( 16)-( 17)

2: μs = max(μ s , μk ). Output μk ← lim ∞ μk (t) μs ← lim ∞ μs (t).
The advantages and drawbacks of this adaptive observer are summarized in Table I.

IV. SECOND APPROACH: RECURSIVE DYNAMICS

FRAMEWORK

The second approach we propose corresponds to a recursive dynamics interconnection framework. Such a framework is particularly suited for sensing and estimation since it allows structuring the model on some simpler and more realistic interconnected dynamical subsystems (blocks) and better exploiting the interconnection structure. The method has been developed in [START_REF] Auriol | Sensing and computational frameworks for improving drill-string dynamics estimation[END_REF] to estimate the nature of the rock interacting with the drill string (only the axial motion was considered). We believe it can be extended to torsional oscillations to estimate the friction parameters. Let us mention that if all the points of the drill string reach the kinematic mode, then the function F(t, x), defined by (4), does not (directly) depend on time anymore. Thus, similar to what has been done in [START_REF] Auriol | Sensing and computational frameworks for improving drill-string dynamics estimation[END_REF], it becomes possible to apply successive backstepping transformations to each section of the drill string to simplify the structure of each block. Denote by x i the spatial coordinate of the junction point between two sections of the drill string. Using the method of characteristics and the continuity of torque and velocity, we can express the state at x = x i as a function of delayed and future state values at x = x i+1 . More precisely, we have

ω i (t, x i+1 ) = f i (ω i (•, x i ), τ i (•, x i )) + h i ω (18) τ i (t, x i+1 ) = g i (ω i (•, x i ), τ i (•, x i )) + h i τ , (19)
where f i and g i are two functions that only depend on the physical parameters of the drilling device (except µ k ) and past and future values of ω i (•, x i ) and τ i (•, x i ), the largest time-shift being fixed by the parameters of the system. The functions h i ω and h i τ are constant functions that linearly depend on µ k . Computations to obtain explicit expressions for these different functions are quite involved and cannot be given here due to space restrictions. They follow the same path as the ones given in [START_REF] Auriol | Sensing and computational frameworks for improving drill-string dynamics estimation[END_REF]. Using equations ( 18) and ( 19), we can express τ (t, L) as a function of top-drive angular velocity and torque measurements. Since the bit is offbottom, we know that τ (t, L) = 0. Thus, it becomes possible to use classical regression techniques (such as Recursive Least Squares) to estimate the kinematic coefficient µ k . Indeed, if we have M measurement points, this coefficient is estimated by minimizing the following cost function

{μ kin } = argmin µ k M i=1 (τ (t i , L)) 2 , (20)
where τ (t i , L) is the estimated version of the torque obtained using equations ( 18) and [START_REF] Kazemi | Across-domains transferability of deep-red in de-noising and compressive sensing recovery of seismic data[END_REF]. Obviously, the proposed estimation of the torque-on-bi only holds if we have ω > ω c all over the drill string (otherwise, non-linear terms appear, which makes the estimation extremely difficult). Thus, we need to apply equations ( 18) and ( 19) to all points of the drill string (in practice, to a sufficiently large number of points) to verify that it is the case. Among all the estimations of the drill-bit source signature, only the ones for which this condition is fulfilled can be used to estimate the kinematic friction term. Note that this estimation procedure is done in near real-time (contrary to the previous adaptive observer which is done in real-time) due to the fact that ( 18) and ( 19) require future values of the measurements.

To verify that the condition ω > ω c is satisfied, we need the knowledge of ω c . However, knowing a decent upperbound could be sufficient and provide an estimation of ω c . More precisely, let us denote ωc such an upper-bound. Applying the previous procedure for ω > ωc all over the drill string, we can estimate µ k , assuming that enough data points Near-real time algorithm Procedure to estimate ωc µs is not estimated Good robustness properties [START_REF] Auriol | Sensing and computational frameworks for improving drill-string dynamics estimation[END_REF] Computationally expensive (τ i ) verify these conditions (to apply the least-squares algorithm). Note that we can modify the reference set-point in the PI control law to guarantee this last condition. Then, since µ k is now known, we can reiterate the estimation procedure while decreasing the bound ωc . Once the quadratic error M i=1 (τ (t i , L)) 2 becomes larger than a chosen threshold , it means that we have reached ω c since the non-linear static friction terms are now modifying our estimation of the torque on the bit. This procedure is described by Algorithm 2. The values of ω c are updated using a dichotomy approach (tolerance η). The threshold should be large enough to encompass measurement noise. The advantages and drawbacks of this estimation procedure are summarized in Table II.

Algorithm 2 Estimation of µ k and ω c with a recursive dynamics framework. Require: τ (t, 0), ω(t, 0), ωc , ω c , , η Compute τi (t, L) with equations ( 18) and ( 19) for all the data points that verify ω(t, x) > ωc , for all t.

2:

Recursive Least Squares algorithm:

{μ kin } = argmin µ k M i=1 (τ (t i , L)) 2 . 3: Update ω k c : ω c +ωc 2 4:
Compute S = M i=1 (τ (t i , L)) 2 for all the data points that verify ω(t, x) > ω k c .

5:

If S > , ω c ← ω k c . Else, ωc ← ω k c .

6:

If ωc -ω c > η go to Step 3. Output μk ← μk , ωc ← ω k c .

V. THIRD APPROACH: MACHINE LEARNING ESTIMATION

The third approach we present in this paper is based on Machine Learning ( [START_REF] Goodfellow | Deep learning[END_REF]). Hence, it does not require any specific knowledge of the system or an accurate model. Note that adaptive neural networks have already been applied in [START_REF] Fu | Neural network-based finitetime command filtering control for switched nonlinear systems with backlash-like hysteresis[END_REF] or [START_REF] Yu | Neuroadaptive finite-time control for nonlinear mimo systems with input constraint[END_REF] for control purposes. The main idea behind this approach is that the measured top-drive torque and angular velocity (or in closed-loop the output of the PI control law) depend on the geometry and the physical parameters of the drilling device (and in particular of the friction factors). Thus, the spectrum of these signals should contain sufficient information to estimate the static and kinematic friction terms. Let us denote y(t) = τ m (t) the measured motor torque (that is related to ω 0 due to [START_REF] Åström | Feedback systems: an introduction for scientists and engineers[END_REF]). Due to the top-drive sensors sampling rate, this signal is not updated continuously. We will denote Y the vector obtained by concatenating the different values of y at each sampling step. We then compute the discrete Fourier transform of Y to obtain the spectrum of the signal y. We recall that the discrete Fourier transform of a signal Y with M components is a complex function that is defined by

Ŷ (k) = M n=1 y(k)e -2πi * (k-1)(n-1) M . ( 21
)
The spectrum of the signal y is characterized by several attributes: the number of peaks, the dominant peak (highest value of the modulus of Ŷ ), the dominant peak after a given cut-off frequency, the corresponding frequencies.

The values of some of these attributes are related to the physical parameter of the systems. Thus, for a set of known physical parameters that characterize the well (including its geometry), we can run thousands of simulations, only modifying the parameters we want to estimate (namely µ k and µ s ) between each simulation. This gives us a set of data (known as training set) for which the correct values of the friction parameters are known. Our algorithm will learn from this data set and find the suitable correlations between the previously defined attributes and the unknown parameters.

Once adequately trained, the machine can be applied to make accurate predictions for new data sets. We choose to develop a neural network algorithm using a sigmoid activation function [START_REF] Goodfellow | Deep learning[END_REF], [START_REF] Nwankpa | Activation functions: Comparison of trends in practice and research for deep learning[END_REF]. The network is trained by simulating thousands of test points for which we know the correct values of the friction parameters. From these simulations, we generate the spectrum of the corresponding output signals and compute several relevant features. While training our neural network, it appeared that the relevant attributes for learning are the two dominant gains (and their corresponding frequencies). We then design our neural network (i.e., we choose the number of neurons, layers, weights) and train it using cross-validation. The algorithm is now ready to be applied to an unknown dataset. The details of the algorithm are summarized in Algorithm 3. The Algorithm 3 Estimation of µ s and µ k using a ML algorithm

Require: τ m (t), training dataset (τ k m (t), µ k k , µ k s ) Training of the network 1:
Compute the Fourier transform of τ k m (t) using (

Compute the two dominant gains and the corresponding frequencies

3:

Using cross-validation, design a neural network with sklearn.neural network.MLPRegressor Estimation

4:

Compute the Fourier transform of τ m (t)

5:

Compute the two dominant gains and the corresponding frequencies

6:

Run the neural network with this new point Output μs ← μs , μk ← μk advantages and drawbacks of this machine learning algorithm are summarized in Table III.

The generalizability of a convolutional neural network depends on the number of learnable free parameters and the depth of the neural network [START_REF] Baum | What size net gives valid generalization?[END_REF], [START_REF] Zhou | Universality of deep convolutional neural networks[END_REF]. In most applications, the comprehensive dataset is not available and the learned network cannot properly model the complexities of the system. To remedy this shortcoming, one needs to take advantage of active learning, model-based learning, and transfer learning to make the learned mapping function general enough so that the learned operator gives proper parameter estimation when applied on a new dataset [START_REF] Cohn | Improving generalization with active learning[END_REF], [START_REF] Yosinski | How transferable are features in deep neural networks?[END_REF], [START_REF] Tzeng | Simultaneous deep transfer across domains and tasks[END_REF], [START_REF] Kazemi | Across-domains transferability of deep-red in de-noising and compressive sensing recovery of seismic data[END_REF].

VI. SIMULATION RESULTS We now test our three estimation algorithms against the simulation model given in [START_REF] Aarsnes | Axial and torsional self-excited vibrations of a distributed drill-string[END_REF] using the wellbore survey from [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF], where the well is vertical for 300m, builds to 40 degrees inclination and holds until to a measure depth of 1250m and then builds and holds at 45 degrees for the remainder of the well.

A. Comparison of the three algorithms on two test cases

We first compare the three estimation algorithms on two test case studies. All the numerical values for the different parameters can be found in [START_REF] Aarsnes | Axial and torsional self-excited vibrations of a distributed drill-string[END_REF]. We consider a reference set-point for the top drive angular velocity of 60 RPM. In the first case study, the side forces static friction term is chosen to be equal to 0.42, while the kinematic friction term is equal to 0.29. In the second one, we have µ s = 0.45 and µ k = 0.28. In both cases, we have ω c = 14 RPM. We assume that the measured signals are subject to a white Gaussian noise that is characterized by its signal-to-noise ratio (SNR) [START_REF] Kazemi | Surface-consistent sparse multichannel blind deconvolution of seismic signals[END_REF]. Seventy-five simulations and estimations have been performed for each algorithm. The SNR for all the signals is equal to 10. We give in Table IV and Table V the mean and standard deviation for the estimations given by each algorithm.

Our three algorithms allow satisfying estimations of the friction parameters. The two first algorithms are modelbased and consequently require a reliable model of the drill string dynamics (which is not the case of the machine learning algorithm). Although being the most accurate, the Recursive dynamics framework-based algorithm cannot (in its current form) provide an estimation of the static friction term. Finally, even if the Machine Learning algorithm shows satisfactory performance, it requires sizeable training datasets, which can be time-consuming to both generate and then subsequently train.

B. Machine Learning Results

To extend the usefulness of the data-driven model, we considered expanding the training dataset to include a range of bit depths, leading to a range in static and kinematic torques. Using two training data sets, we evaluate the performance of the neural network on a test data set generated for the considered well. Training data set 1 (DS1) contains 1000 data points with parameters ranges of µ s ∈ [0.2, 0.83], f rat ∈ [0.2, 0.65] (where f rat is such that µ k = f rat µ s ) and for well depths ranging from 2,500 m to 2,770 m. Training data set 2 (DS2) contains 1000 data points with the same ranges for µ s and f rat , however, well depth ranges from 1,700 m to 2,870 m. The test data set (TS) contains 125 data points with µ s ∈ [0.210.76], f rat ∈ [0.23, 0.6] and measure depths ranging from 1,760 m to 2,752 m. Using the algorithm described in Algorithm 3, a neural network is trained on DS1 (NN1) and on DS2 (NN2) and then applied to TS.

As the test data set contains depths outside of the range of DS1, the performance of NN1 is poor, as expected. Bit depth and the weighted average of the inclination of the drill pipe are used as additional input parameters. However, the neural network cannot generalize the predictions as no implicit knowledge of drill string dynamics is embedded in the neural network. Figure 2 shows the actual vs predicted friction values for the training data set, with static friction shown in blue and kinematic friction shown in red.

NN2 performs well when applied to the data in TS as the TS2 spans the input parameter range, for µ s , µ k and bit depth. It is noted that estimates of static friction are worse than the estimated of kinematic friction, as is expected. Figure 3 shows the actual vs predicted values of friction for the training data set.

This shows that the machine learning model is generalizable, so long as the training data set spans a sufficient parameter space. Future work in embedding some of the basic physical properties of the drill string system into the neural network may improve performance outside of the training parameter space. Some preliminary work has been in this space for seismic data processing and inversion shows promise for application to problems governed by the wave equation [START_REF] Sun | Physics-guided deep learning for seismic inversion with hybrid training and uncertainty analysis[END_REF].

VII. CONCLUDING REMARKS

In this paper, we have designed three algorithms that provide real-time or near real-time the friction corresponding to the side-forces on a drilling device. Estimating these parameters is a well-known challenge in the drilling industry and, if understood, presents the opportunity to improve drilling operations through feedforward or model predictive control. It is a crucial step to develop stick-slip mitigation strategies. The three proposed procedures correspond to an adaptive observer (originally introduced in [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF]), to an algorithm based on recursive dynamics framework, and to a machine learning algorithm. All of them only require surface measurements, making them realistically implementable on drilling devices. Their specific advantages and drawbacks have been emphasized. These algorithms have then been implemented on a previously validated torsional drill string model, with specific attention paid to the machine learning procedure. This algorithm is easy to implement and fast to run once properly trained but may not be robust to changes in the system parameters (as a new long training phase would be necessary). The three algorithms provide satisfying estimations of the friction parameters. Future works should include testing against field data and an analysis of the computational effort for each procedure. The expected outcome is a hybrid strategy to estimate friction parameters in drilling.

Fig. 1 .

 1 Fig. 1. Schematic indicating the distributed drill string of length L lying in deviate borehole. Inset: Schematic illustrating the friction source term S(ω, x). The shaded region represents the angular velocities for which a constant value of static torque is assumed and the red curve indicates the dynamic torque as a function of angular velocity.

Fig. 3 .

 3 Fig. 3. Comparing real friction versus predicted friction for NN2, where the training data set includes the range of depths included in the test data set. Performance is significantly improved, particularly for kinematic friction.

TABLE IV MEAN

 IV AND STANDARD DEVIATION (IN PARENTHESIS) OF THE FRICTION FACTORS ESTIMATION USING THE THREE PROPOSED ALGORITHMS FOR THE FIRST TEST CASE (µs = 0.34 AND µ k = 0.187).

	Algorithm	μs	μk
	Adaptive observer	0.358 (0.02)	0.192 (0.01)
	Recursive dynamics framework	N.A	0.1892 (0.03)
	Machine learning algorithm	0.332 (0.023)	0.19 (0.01)

TABLE V MEAN

 V AND STANDARD DEVIATION (IN PARENTHESIS) OF THE FRICTION FACTORS ESTIMATION USING THE THREE PROPOSED ALGORITHMS FOR THE SECOND TEST CASE (µs = 0.45 AND µ k = 0.28).Fig. 2. Comparing real friction versus predicted friction for NN1, where the training data set has a limited range in depths. Poor performance is observed as the neural network fails to properly predict friction for depths outside of the training data set range.

			1
			0.8
			0.6
			0.4
			0.2
			0
			0	0.2	0.4	0.6	0.8	1
	Algorithm	μs	μk
	Adaptive observer	0.50 (0.015)	0.289 (0.01)
	Recursive dynamics framework	N.A	0.287 (0.03)
	Machine learning algorithm	0.463 (0.03)	0.285 (0.02)


Benjamin Mathon
email: benjamin.mathon@gipsa-lab.inpg.fr

Patrick Bas
email: patrick.bas@gipsa-lab.inpg.fr

François Cayre
email: francois.cayre@gipsa-lab.inpg.fr

Fernando Pérez-González

Distortion Optimization of Model-Based Secure Embedding Schemes for Data-Hiding

Keywords: Watermarking, Security, Hungarian method, Distortion Optimization

This paper is the continuation of works about analysis of secure watermarking schemes in the case of WOA (Watermarked Only Attack) framework. In previous works, two new BPSK spread-spectrum watermarking modulations, Natural Watermarking (NW) and Circular Watermarking (CW), have been proposed and have been shown to be more secure than classical modulations. Because security is guaranted using specific distributions of watermarked contents, we propose to use optimal model-based embedding to ensure security while minimizing the overall distortion. Additionally, we propose a new secure watermarking scheme based on distribution of vector norms in the Gaussian case. We illustrate model-based embedding performance in the case of Gaussian signals and show that this approach not only allows to achieve excellent level of security in the WOA framework, but also allows to minimize distortion. Finally, a comparison of the robustness of the proposed embedding schemes is performed.

Introduction

Watermarking is a mean to hide information into digital contents (images, sounds, videos). This hidden information can be used for copyright or copy protection applications, integrity checking, or fingerprinting in order to control each copy of a numerical document. Our works focus on copyright and copy protection. The embedding of the message must meet many constraints:

imperceptibility: the hidden information must not impair the original content for regular use, capacity: in multi-bit watermarking, one must ensure a sufficient number of bits can be reliably hidden into the host content, robustness: the hidden message should still be readable after common media processing, security: this last constraint can be viewed as "the inability by unauthorized users to access, remove, read or write the hidden message" [START_REF] Kalker | Considerations on watermarking security[END_REF]. Security in general is based on Kerckhoffs' principle [START_REF] Kerckhoffs | La cryptographie militaire[END_REF]: the key is the only unknown parameter for the adversary.

This secret key is used to embed and detect (zero-bit watermarking) or decode (multi-bit watermarking) the watermark. Security is different from robustness. Robustness attacks relate to watermark survivability under common processing (in the case of still images, one may want to resist geometrical deformations, compression or noise addition). These attacks are generally not intentional.

On the contrary, security attacks are intentional and relate to the estimation of a part or all the secret key [START_REF] Comesaña | Fundamentals of data hiding security and their application to spread-spectrum analysis[END_REF]. In multi-bit watermarking, the key is defined by the location of a set of codewords in a subspace. To embed a message in a host content, it must be placed in the decoding region of the right codeword. If an adversary learns the secret key, he can alter the message while minimizing the attack distortion with a 100% probability: his attack becomes deterministic. For example, he can design an attack in order to "push" the watermarked content into the nearest (wrong) decoding region, he can read the hidden message and copy the watermark in another content (copy attack [START_REF] Kutter | The watermark copy attack[END_REF]).

This paper deals with the case of WOA (Watermarked Only Attack) [START_REF] Cayre | Watermarking security: Theory and practice[END_REF], the adversary has only access to several marked contents and to the source code of the watermarking algorithm (Kerckhoffs' principle). One should notice that the adversary can model the distribution of the host contents (Gaussian mixture or generalized Gaussian distribution for DCT or wavelet coefficients for example) since he knows the watermarking space. The adversary's goal is twofold. On one hand, he wants to model the conditional distribution of marked contents given the secret key, and on the other hand, he wants to estimate each codeword location.

Based on these ideas, we can find [START_REF] Cayre | Kerckhoffs based embedding security classes[END_REF] the definition of four security classes to rank watermarking schemes security in the WOA framework (see Fig. 1 for relationships among them).

The first class is insecurity. In this class, the conditional distribution (given the key) of marked contents is not the same for all keys. By exhaustive search (or a more involved technique), he can estimate both the private subspace and the codewords.

A watermarking scheme belongs to the second class, key-security, if, for a subset of keys, the conditional distribution of the marked contents given each key is the same. The adversary can find this subset of keys and he can find the secret subspace but he cannot gain more information about the codewords.

The third class is called subspace-security. In this case, the conditional distribution of marked contents given the key will be the same for all keys. Therefore the pirate cannot gain any information about the key (he has not access to the secret subspace).

The last class stego-security relates to steganography: the distribution of marked contents is the same than that of the host contents. The adversary cannot decide whether the contents are marked or not. Obviously, security in the WOA framework is strongly linked with the distribution of the contents after watermark embedding. In [START_REF] Mathon | Practical performance analysis of secure modulations for woa spread-spectrum based image watermarking[END_REF], authors have applied two secure modulations (namely Circular and Natural Watermarking, resp. CW and NW in the sequel) to still images. These modulations are used to modify the distribution of marked contents in order to be stego-secure or subspace-secure for NW and key-secure for CW. Note that such an approach is similar to the one proposed by Sallee in steganography [START_REF] Sallee | Model-based steganography[END_REF], where perfect secrecy [START_REF] Shannon | Communication theory of secrecy systems[END_REF] is guaranteed by constraining the distributions of stego contents to be identical with the distributions of cover contents.

In this paper, we propose a new method to watermark signals in order to fit a chosen target distribution in an optimal way. This method uses the Hungarian algorithm, which minimizes distortion on average (in the sense of Euclidean distance) between points of host distribution and points of the target distribution. Section 2 recalls basics on BPSK-based SS watermarking schemes with an emphasis on unsecure and secure modulations and presents a watermarking scheme based on distribution of norm of signals. Section 3 presents our model-based embedding scheme by using Hungarian method. Section 4 compares our implementations of classical versus model-based embedding from security, distortion and robustness point of view over 2000 Gaussian signals.

Secure watermarking schemes

This section recalls definitions of unsecure and secure modulations, moreover we propose also a new secure embedding scheme based on the modification of the norm of the host vector.

Notations and definitions

We first list the conventions used in this paper. Data are written in small letters. Vectors and matrices are set in bold fonts. Vectors are written in small letters and matrices in capital ones. x(i) is the i-th component of a vector x and x j is the vector x associated to a j-th observation. We write (x(0), x(1), x(2), ...) the content of a vector x. We note [a;b], ]a;b[, [a;b[ and ]a;b] real-intervals. Sets are noted in capital letters and vect (A) represents the vector space generated by A. p (x j ) denotes the distribution of vectors x j . σ 2

x denotes the variance of a signal x and .|. denotes the usual scalar product.

We want to hide a message m of N c bits in a host vector x ∈ R Nv . So we create a watermark signal w ∈ R Nv in order to obtain y = x+w the watermarked signal. The secret key K used to embed and decode the message m is the seed of a PRNG. With K, we generate N c Gaussian carriers u i ∈ R Nv . Each carrier is able to hide one bit. In order to have a null ISI (Inter Symbol Interference), carriers must be orthogonal. Thanks to s : {0, 1} → R, a modulation, we can create w by:

w = Nc-1 i=0 u i s(m(i)). (1)
Distortion is assessed by means of the W CR (Watermark-to-Content Ratio):

W CR [dB] = 10 log 10 σ 2 w σ 2 x . (2)
We model robustness attacks by adding Gaussian noise n. So we consider the attacked vector r = y + n. Attack strength is assessed by means of the W CN R (Watermarked Content-to-Noise Ratio):

W CN R [dB] = 10 log 10 σ 2 y σ 2 n . (3)
Decoding is classically obtained by correlations z:

z r,ui = Nv-1 j=0 r(j)u i (j). ( 4
)
We consider m the estimated message, so we have for each bit:

m(i) = 0 if z r,ui > 0, 1 if z r,ui < 0. ( 5
)
We measure robustness of the watermarking scheme by BER (Bit Error Rate) between the estimated and the original message:

BER(m, m) = 1 N c Nc-1 i=0 m(i) ⊕ m(i). (6)
For a pirate, there is no difference between estimating the carriers or getting K in the WOA context. According to [START_REF] Cayre | Kerckhoffs based embedding security classes[END_REF], the security of a watermarking scheme relies on the properties of the conditional distribution p(y j |K). In the case of spead-spectrum techniques, points of conditional distribution given the carriers are the N c -tuples (z y,u0 , ..., z y,u Nc -1 ) in the subspace spanned by the carriers.

Unsecure SS modulations:

Classical modulation SS (Spread Spectrum) is given by:

s SS (m(i)) = γ(-1) m(i) . ( 7
)
This modulation is analog to BPSK modulation. Parameter γ is used to set the power of watermark. It is a function of W CR. ISS (Improved Spread Spectrum) [START_REF] Malvar | Improved spread spectrum: a new modulation technique for robust watermarking[END_REF] uses side-information to improve robustness:

s ISS (m(i)) = α(-1) m(i) -λ x|u i u i 2 , (8)
where .|. denotes the usual scalar product, α and λ are computed to achieve host-interference rejection and error probability minimisation given Noise-to-Content power Ratio:

N CR [dB] = 10 log 10 σ 2 n σ 2 x , (9)
where n denotes Gaussian noise. Previous works [START_REF] Cayre | Watermarking security: Theory and practice[END_REF] have shown that SS and ISS are unsecure, and that carriers estimation is possible (see tests on [START_REF] Mathon | Practical performance analysis of secure modulations for woa spread-spectrum based image watermarking[END_REF]).

Secure modulations:

NW (Natural Watermarking) [START_REF] Cayre | Kerckhoffs based embedding security classes[END_REF] modulation is defined by:

s N W (m(i)) = (-1) m(i) x|u i | x|u i | -1 x|u i u i 2 . ( 10
)
NW belongs to the so-called stego-secure class and is suitable for steganography applications. CW (Circular Watermarking) [START_REF] Cayre | Kerckhoffs based embedding security classes[END_REF] modulation is defined by:

s CW (m(i)) = α(-1) m(i) d(i) -λ x|u i u i 2 , (11)
where α and λ are computed the same way than with ISS and d is generated at each embedding as follows from g ∼ N (0, 1), this parameter is used to randomly spread the correlations of the mixed signals on the whole decoding regions:

d(i) = |g(i)| g . ( 12
)
This parameter d enables the following property of circularity :

p z y,u0 , ..., z y,uN c -1 = p   Nc-1 i=0 z 2 y,ui   . ( 13
)
The circularity of the distribution allows us to say that for a subset of several keys (all bases of vect({u i })), the distribution of marked signals will be the same. CW belongs to the so-called key-secure security class.

A new secure embedding based on the χ 2 distribution

Based on previous security assessment, we are able to propose a stego-secure watermarking scheme χ 2 W (CHI2 Watermarking) by modifying norm of Gaussian host signals after embedding, while keeping the same distribution between original and watermarked contents. Distribution of these norms can be modeled by a χ 2 law. So we define codewords location by real-intervals in the set of norms of Gaussian vectors. Finally, in order to embed a secret message, we chose randomly a norm in the corresponding interval and we multiply the host vector in order to have the desired norm. Contrary to BPSK modulations, the watermarking subspace (space of norms) is not private, secret relies only on the partition of the real-positive axis representing the norms. There is no security flaw because we works on the WOA framework. Adversaries do not know the embedded message. However, this embedding is easy to implement, it enables to achieve stego-secure embedding since the distributions of original and watermarked contents are the same, and brings another scheme to compare with in term of robustness. We use previous notations of Part. 2.1. We want to create a watermarked signal y = αx, α ∈ R + . This method is based on the distribution of the norms of host Gaussian vectors. In fact, if x ∈ R Nv with x ∼ N (0, 1),

x 2 ∼ χ 2 (N v ) (Chi-2 law of degrees N v ), .
representing the euclidean norm. Codewords are sets in a partition of [0, +∞[. To embed a secret message m in a host vector x, we randomly choose a norm y 2 in the corresponding real-interval and we compute:

y = y 2 x 2 x. (14)
We obtain the watermark signal w = yx. This process can be considered as a variant of Moulin and Briassouli stochastic embedding [START_REF] Moulin | A stochastic qim algorithm for robust, undetectable image watermarking[END_REF] who work on different host distributions. Different consequences arise from this new embedding scheme:

-The means of choosing a norm in the right codeword is not optimal (from the distortion point of view): we generate real numbers until we have one of them in the desired interval. -To define the secret partition on the real-positive axis, we use an estimator of the quantile function of χ 2 distribution. -With the condition of message equiprobability, we must have more than one codeword for each message. Without this condition, a pirate can find the secret partition by using a quantile function (he separates the p.m.f. into 2 Nc parts of probability 1 2 Nc ). We denote N w the number of codewords used. Fig. 2 shows a secure partition of real positive axis and the associated distribution with parameters Note that the decoding regions are the same for NW and CW (they are delimited by hyperplanes) but different from the ones related to χ 2 W (delimited by hyperspheres). Fig. 3 shows 2D representations of coding regions for spread-spctrum schemes and for χ 2 W. We use N c = 2.

N c = 2, N w = 8, N v = 55.

Minimisation of the embedding distortion

We have seen that the level of security of the previously watermarking methods is given by the distribution of the signals after watermarking. Note that these distributions are defined by the distribution of correlations for BPSK modulations and the distribution of the norms for χ 2 W. However, from the distortion point of view, these implementations are not optimal. We propose in the next section a new scheme which can ensure a given distribution of our marked signals while minimizing the embedding distortion. We want to associate each point of host distribution with each point of a chosen distribution with a minimal average Euclidean distance between these two points (distortion is proportional to the distance). The Hungarian algorithm is a mean to solve this problem. To explain this algorithm, some reminders about graph theory are useful.

Minimal cost perfect matching in a bipartite graph

A bipartite graph is a graph G = (V, E) with the following property: there exists a partition E,P ) is a bipartite graph where each edge is weighted by a function P : E → R. A perfect matching M of G is defined as a subset of E with N m elements where each vertex is incident with exactly one member of M . In this paper, we are interested in the Assignement Problem (AP), we search the minimal cost perfect matching M * , i.e., a perfect matching whose the sum of weights of edges is minimal. More precisely we search:

V = A ⊔ B, each edge of E is of the form [a, b] with a ∈ A and b ∈ B. Moreover G satisfy |A| = |B| = N m . A weighted bipartite graph G = (V,
M * = arg min M t∈M P (t). (15)
Fig. 4 shows a weighted bipartite graph and its minimal cost perfect matching.

The Hungarian method for the AP

The Hungarian method [START_REF] Kuhn | The Hungarian method of solving the assignment problem[END_REF] is an efficient algorithm to solve the AP in a weighted bipartite graph in polynomial time (O(N m 3 )). We consider G = (V, E, P ) a weighted bipartite graph with:

-V = A ⊔ B, -A = {a 0 , ..., a Nm-1 }, -B = {b 0 , ..., b Nm-1 }.
We consider D ∈ M Nm,Nm (R), a matrix initialized with D(i, j) = P ([a i , b j ]). The goal is to choose N m elements of this matrix in order to have each row and each column containing one chosen element. In fact, the minimal cost perfect matching is the set of edges corresponding to these choosen elements. The Hungarian algorithm does the following:

1. Subtract the entries of each row by the row minimum: each row has at least one zero, all entries are positive or zero. 2. Subtract the entries of each column by the column minimum: each row and each column has at least one zero. 3. Select rows and columns across which to draw lines, in such a way that all the zeros are covered and that no more lines have been drawn than necessary.

A test for optimality:

-If the number of the lines is n, choose a combination from the modified cost matrix in such a way that the sum is zero. -If the number of the lines is < n, go to 5. 5. Find the smallest element which is not covered by any of the lines. Then subtract it from each entry which is not covered by the lines and add it to each entry which is covered by a vertical and a horizontal line. Go back to 3.

So we have:

M * = {[a i , b j ] : D(i, j) = 0}. ( 16
)

Application for NW and CW embedding

Construction of bipartite graphs: We want to create 2 Nc bipartite graphs which contain, for the host partition, points of distribution of several host signals (correlations). We construct points of target distributions by selecting only points in codeword of the desired message. The goal is to find the minimal cost (euclidean distance) perfect matching between the two partitions. Formally, we use N m host signals {x i } i=0,...,Nm-1 . For each host signal x i we construct the correlation host vector:

z xi = z xi,u0 , ..., z xi,uN c -1 .
We want to construct 2 Nc target distributions. As we have seen before, NW and CW modulations can permit to construct distributions we want to obtain in order to match the security class we want (stego-secure, key-secure, ...). We proceed by watermarking for each message m, host signals x i with the chosen modulation to obtain y i . Finally, we construct the correlation marked vector z yi = z yi,u0 , ..., z yi,uN c -1 . We obtain, for each message, the weighted bipartite graph G = {X ⊔ Y, A, P } with:

-X = {z xi } i=0,...,Nm-1 , -Y = {z yj } j=0,...,Nm-1 (depends on m), -A defines the set of edges [z xi , z yj ] of the graph G, -P is the weight function of the edges of G. P ([z xi , z yj ]) = z xi -z yj 2 .
We obtain 2 Nc minimal cost perfect matchings M * k between host correlations and marked correlations by using Hungarian method.

Mapping reduction: In previous section, we construct one bipartite graph for one message to embed. In order to reduce complexity of bijections calculus, we can use property of symmetry of our distributions (the axis of symmetry are the carriers); points of target distributions are computed in order to embed a constant message (for example the message (1, 1, ..., 1)). For the rest of this article, we notate N m -map a triplet (X, Y, M * ) constructed with N m host signals. Fig. 5 shows a correlation bipartite graph with N m = 3, N c = 2 and its minimal cost perfect matching found by Hungarian method. Model-based embedding: We consider (X, Y, M * ) the N m -map constructed in the previous section. Now, we want to mark a signal x with any message m using the map. First, we compute the correlation host vector z x = z x,u0 , ..., z x,uN c -1 .

We want to associate z x with a point of Y by using M * . We have seen that elements of Y have been constructed in order to embed the message (1, ..., 1). We note indices of m where the bit is different from 1. Note that a matching from an original content to a given codeword enables to generate matchings to any codewords by symmetrising both points along appropriated axes. Consequently sign changes must be made on the coefficients of z x in the indices that have undergone symmetries. Afterwards inverse symmetry must be performed after watermarking in order to embed the correct message m. Formally, we construct Rz x ,

R ∈ M Nc,Nc (R), R(i, j) = 0 if i = j, (-1) m(i)+1 if i = j.
Next, we find the nearest neighboor (minimal euclidean distance) of Rz x in X, for example z xi 0 . Thanks to the perfect matching M * we find z yj 0 , the correspondance of z xi 0 . Next, we apply inverse symmetry to compute the correlation marked vector z y :

z y = R -1 z yj 0 = Rz yj 0 .
So, we obtain the correlation vector of our marked signal. By a difference between z y and z x , we have z w , the watermark correlation vector. Proper retro-projection of this signal in the N v D-space is assured by :

w = Nc-1 i=0 z w (i) u i |u i u i .
And finally, we compute the watermarked signal y = x + w. Fig. 6 shows this process by using the 3-map constructed, see Fig. 5.

Application for χ 2 Watermarking

Construction of bipartite graphs: We want to create 2 Nc bipartite graphs which contain, for the host distribution, norms of several host signals. We construct points of target distribution by selecting only real points in codewords of the desired message. We generate N m Gaussian vectors x i and, for each message, N m Gaussian vectors y j with y j 2 in the right codewords. We can construct, for each k = 0, ..., 2 Nc -1 a bipartite graph G = {X ⊔ Y, A, P } with:

-X = { x i 2 } i=0,...,Nm-1 , -Y = { y j 2 } j=0,...,Nm-1 , -A defines the set of edges [ x i 2 , y j 2 ]
of the graph G, -P is the weight function of the edges of G. We use the absolute value of the differences of norms. By using the Hungarian algorithm, we find 2 Nc minimal cost perfect matching M * k between host norms and marked norms (functions of embedded message). For the rest of this article, we denote N m -smap the set {(X, Y k , M * k )} k=0,...,2 Nc -1 constructed with N m signals.

Embedding: Process is similar to SS model-based embedding. To embed a vector x with message m and a N m -smap; we calculate x 2 . Next we find the nearest neighboor of x 2 in X. By M * k (k depends on m), we find y 2 . Finally, we obtain y = y 2

x 2 x.

Experiments

The goal of this section is to assess the preservation of the distributions after the Hungarian method, the impact of this method on distortion and the general robustness of the three secure embedding schemes we presented.

Numerical values and assessments

In practice, N c = 2. For NW and CW, we use N v = 256 and for χ 2 W, we use N v = 55 (in order to have the same distortion for the three schemes, W CR = -18dB), N w = 8. Tests are made with 2000 host Gaussian signals. For χ 2 W, in order to have the equiprobable condition, we use an estimator of a fractile function of χ 2 distribution given in [START_REF] Hill | Algorithm 299[END_REF] which uses an estimator of the repartition function of a normal distribution in [START_REF] Ibbetson | Algorithm 209[END_REF], we use the partition of the real-positive axis defined on Tab. 1. We have constructed a 10000-map and a 10000-smap with our signals database by using CW, NW and χ 2 W and we have marked the initial 2000 signals by using on the one hand the secure modulation, and on the other hand the corresponding model-based found by Hungarian method (HCW, HNW and Hχ 2 W). robustness of χ 2 W can be given by the fact that for this embedding, the decoding regions are always very close to each other (see Fig. 3). Consequently one watermarked vector corrupted by noise will have a higher probability to change of coding regions for χ 2 W than for NW or CW. Note that we cannot show robustness of the six schemes with the same distortion because NW modulation does not allow to set a target distortion.

Conclusion

The goals of this paper are twofold: to propose and compare secure embedding schemes for data-hiding and to propose a general method to minimise the global embedding distortion for each scheme. The first point is addressed by proposing the χ 2 embedding scheme which is more fragile than other scheme like NW or Table 2. Distortion for NW, CW and χ 2 W on initial embedding schemes and after using the Hungarian optimisation scheme.

CW. The optimisation regarding robustness appears not to be straightforward and future works will be devoted to find more robust schemes and to find coding regions that will improve the robustness of χ 2 W. Moreover, we found that the Hungarian method is the ideal practical tool to minimize distortion while guaranteeing a given class of security.

Acknowledgments

Fig. 1 .

 1 Fig. 1. Security classes in WOA framework.

Fig. 2 .

 2 Fig. 2. Secret partition of real-positive axis and associated distribution with parameters Nc = 2, Nw = 8, Nv = 55, each message is coded into 2 codewords. We have constructed the bins by generating 2000 Gaussian vectors for each message and calculating their respectives norms.

Fig. 3 .

 3 Fig. 3. 2D representations of coding regions for NW and CW (left part) and for χ 2 W (right part) for 2 bits.

Fig. 4 .

 4 Fig.[START_REF] Kutter | The watermark copy attack[END_REF]. Example of a weighted bipartite graph with two partitions of three vertices and minimal cost perfect matching found by Hungarian Method (weights between two vertices are noted on each corresponding edge).

Fig. 5 .

 5 Fig.5. Projection over two carriers (Nc = 2): correlation bipartite graph construction and minimal cost perfect matching found by Hungarian method with two partitions of three vertices, host correlations and marked correlations (Nm = 3). Euclidean distances (weights) between two vertices are noted on each corresponding edge. The minimal cost perfect matching associates elements of each vertex partition while minimizing the summation of the distances.

Fig. 6 .

 6 Fig.6. Model-based watermarking scheme: illustration of Sec. 3.3. After calculating the correlation vector zx of a host signal, we compute the watermarked correlations zy by using the 3-map of Fig.5(Nm = 3) with the constant message m = (0, 1) (Nc = 2).

Fig. 7 .Fig. 8 .

 78 Fig. 7. (a): distribution of the projections of the host signals over two carriers. (b) and (c): distributions of the projection of the marked signals over two carriers for NW and HNW.

Fig. 9 .

 9 Fig. 9. (a): projection of host signals over the two first components. (b) and (c): projections of marked signals distributions over the two first components for χ 2 W and Hχ 2 W. WCR (classical) WCR (with the Hungarian method) NW -18.07 -20.76 CW -17.97 -19.11 χ 2 W -18.02 -21.65

  Mathon, François Cayre and Patrick Bas are supported (in part) by the European Commission through the IST Programme under Contract IST-2002-507932 ECRYPT and the National French projects Nebbiano ANR-06-SETIN-009, RIAM Estivale and ARA TSAR.

Fig. 10 .

 10 Fig. 10. BER vs WCNR for NW, CW and χ 2 W. For NW, CW and χ 2 W, we have W CR = -18dB. For HNW, HCW and Hχ 2 W, we obtain respectively W CR = -20.76dB, -19.11dB and -21.65dB.

Distribution preservation after the Hungarian method

Fig. 7 shows host, NW and HNW distributions on two carriers. As we can see, distribution of correlations after Natural Watermarking is the same than distribution of host correlations. It is consistent with the definition of stego-security. Moreover, correlations are the same after using Hungarian Method. Our modelbased doesn't impair security. Distribution of CW and HCW on two carriers is shown on Fig. 8, we can see that the distribution is circular and we can conclude that for all bases ( u 0 , u 1 ) of vect(u 0 , u 1 ) the distribution p(y 0 , ..., y 1999 | u 0 , u 1 ) will be the same (rotations of the secret subspace). It is consistent with the definition of key-security, the pirate can access the subspace of the codewords but has no information about the decoding regions. As HNW, HCW doesn't impair security. Fig. 9 shows projections of host, χ 2 W and Hχ 2 W signals over the two first components. The distribution do not change with the two methods. As NW, χ 2 W is stego-secure.

Distortion minimisation

Tab. 2 shows the impact on distortion obtained on average on our 2000 signals for NW, CW, χ 2 W, HNW, HCW and Hχ 2 W. We can see that we gain 2.7dB of distortion for NW, 1.1dB for CW and 3.6dB for χ 2 W. This last result is due to the fact that there are two codewords for on message and these codewords are away in the real-positive axis.

Robustness

Beside distortion and security, the last constraint to assess is the general robustness of the presented schemes. We measure robustness of these secure modulations with and without distortion optimisation. Fig. 10 shows Bit Error Rate functions of chosen W CN R over 2000 signals and we can verify that distortion optimisation does not modify robustness of our schemes. As we can see, CW is more robust than NW which is more robust than χ 2 W. An insight of the poor