Delphine Sinoquet
email: delphine.sinoquet@ifp.fr

Gregory Rousseau

Yohan Milhau

Design optimization and optimal control for hybrid vehicles

Keywords: Hybrid vehicles, Optimal control, Dynamic programming, Pontryagin's principle

In the context of growing environmental concerns, hybrid-electric vehicles appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions.

This paper presents a parametric study focused on variations of the size of the powertrain components, and optimization of the power split between the engine and electric motor with respect to fuel consumption. To take into account the ability of the engine to be turned off, and the energy consumed to start the engine, we consider a second state to represent the engine: this state permits to obtain a more realistic engine model than it is usually done. Results are obtained for a prescribed vehicle cycle thanks to a dynamic programming algorithm based on a reduced model, and furnish the optimal power repartition at each time step regarding fuel consumption under constraints on the battery state of charge, and may then be used to determine the best components of a given powertrain.

To control the energy sources in real driving conditions, when the future is unknown, a real-time control strategy is used: the Equivalent Consumption Minimization Strategy (ECMS). In this strategy, the battery is being considered as an auxiliary reversible fuel reservoir, using a scaling parameter which can be deduced from dynamic programming results. Offline optimization results and ECMS are compared for a realistic hybrid vehicle application.

Introduction

Growing environmental concerns coupled with concerns about global crude oil supplies stimulate research on new vehicle technologies. Hybrid-electric vehicles (HEV) appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions [START_REF] German | Hybrid powered vehicles[END_REF]. Indeed, it is well known that the kinetic energy gained by the classical vehicle after an acceleration can not be recuperated when braking. One of the advantage of the hybrid vehicle is its ability to recuperate a part (depending of efficiencies) of this energy, to store it in a battery, and to re-use it latter. This electric energy can therefore be used by an electric motor to power the vehicle and to save fuel.

Among all the possible configurations for typical HEV, the parallel powertrain is one of the most effective. In a parallel HEV, an electric machine furnishes a traction power alongside with the engine. The benefits of such system come from its ability to provide a traction force to the drive train either with the engine, either with the only electric motor, or with the both. It is thus possible to power the vehicle only with the electric motor, while the engine is being off. Such powertrain configuration has often been studied to minimize fuel consumption (see [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF][START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF][START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF][START_REF] Paganelli | Equivalent consumption minimization strategy[END_REF]. Nevertheless, extra fuel consumption at engine starts is never taken into account; we propose in the later to define a new state to deal with this start-up cost.

The control of hybrid powertrains is more complex than control of classic engine. Actually, the control laws have to deal with the state of charge of the battery, which provides the level of remaining electric energy, and with the variable efficiency of each element of the powertrain. To find the global optimal solution, control techniques such as linear programming, optimal control [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF]Guzzella and Sciarretta 2005), and especially dynamic programming (DP) [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF][START_REF] Guilbaud | Control optimization for a power-split hybrid vehicle[END_REF]Liu and Peng 2006;[START_REF] Sundström | Optimal hybridization in two parallel hybrid electric vehicles using dynamic programming[END_REF]) have been studied. Optimization of energy management strategies on given driving cycles is often used. In general, these techniques do not offer an on-line solution for real-time control, because they assume that the future driving cycle is completely known. Nevertheless, their results can be used as a benchmark for the performance of other strategies, to perform a parametric study, or to derive sub-optimal control laws to be implemented on a vehicle.

The objective of the parametric study is to determine the influence of component sizes on the fuel consumption, and to find the optimal size of some predominant elements for a prescribed driving cycle. In [START_REF] Bonnans | Parametric optimization of hybrid car engines[END_REF], the author considers a bilevel optimization problem, with a simultaneous optimization of design parameters and trajectories. Although the implementation of this method is complicated, and the computation time is high, it allows obtaining the exact optimal size of several elements along a driving cycle. Some other studies [START_REF] Sharer | Vehicle simulation results for plug-in HEV battery requirements[END_REF][START_REF] Assanis | Optimization approach to hybrid electric propulsion system design[END_REF]) consider more accurate models, but do not control engine and electric motor with an optimal controller. In this context, we chose to test different component size values, and to get the best size obtained on an optimal trajectory, thanks to a dynamic programming algorithm.

Among sub-optimal control laws used as real-time controller, ECMS [START_REF] Paganelli | Equivalent consumption minimization strategy[END_REF][START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF]) is based on physical considerations, and is then applicable much more generally than rule-based controllers. This controller is inspired from the Pontryagin Principle, and thus some links can be done between parameters of this controller and dynamic programming results [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]. The performance of this controller resides in the control of the energy ratio (called p in the following) between electrical energy and fuel energy.

In a first part of this paper, we focus on a Citroën C1, equipped with a gasoline engine displacement of 1000 cc. From this vehicle, we consider an additional source of power given by an electric motor located near the wheels. This motor takes its energy from an additional battery. A model of this powertrain has been implemented in AMESim environment for validation purposes, and also in Matlab, for optimization. The later model is a simplified static model that allows fast computations, and includes the start-up cost of the engine. Two different optimization algorithms are presented and applied on this model, depending on the needs:

a Dynamic Programming algorithm [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF][START_REF] Scordia | Systematic elaboration of online energy management laws for hybrid vehicles[END_REF][START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], which provides the optimal control strategy on a prescribed cycle, and the potential of a powertrain. In the algorithm is considered the extra fuel consumption (called L 2 in the following), a sub-optimal control law, based on the Principle of Pontryagin, used as a real-time control law.

2 System modelling and optimal control problem

Optimal control problem

The optimal control problem under study consists in minimizing the fuel consumption of the vehicle along a prescribed vehicle cycle, taking into account physical constraints from battery, engine and electric motor. Let x 1 be a state, corresponding to the state of charge of the battery and x 2 ∈ {0, 1} be a discrete state representing the state of the engine: x 2 = 0 means the engine is off, while x 2 = 1 means the engine is on. We denote by x = {x 1 , x 2 } the state vector.

We define the control u = {u 1 , u 2 , u 3 } ∈ U(x 2 , t) by u = ⎧ ⎨ ⎩ u 1 the power split ratio between the engine and the electric motor, u 2 the gear number, u 3 starting up/stopping the engine,

(1)

with U(x 2 , t) is such that U(x 2 , t) = {(u 1 , u 2 , u 3 )\u 1 ∈ U 1 (x 2 , t), u 2 ∈ U 2 (t), u 3 ∈ U 3 (x 2 )}, ( 2
)
where

U 1 (x 2 , t) = [u min 1 , u max 1 ], U 2 (t)
are discrete values, and U 3 (x 2 ) are discrete values {-1, 0} or {0, 1} depending on x 2 value (see Sect. 3). The control u 1 corresponds to the classical control for HEV, u 2 may be optimized too, especially for automatic transmissions or CVT transmissions (continuous variation of the gear ratio), and u 3 is only linked to x 2 state. Let g(x 1 (T ), T ) be a penalization term handling a constraint on the final state of charge of the battery, and L(x(t), u(t), t) be the instantaneous fuel consumption, given by

L(x(t), u(t), t) = ⎧ ⎨ ⎩ L 1 (u 1 (t), t) if x 2 = 1 (engine started), 0 i fx 2 = 0 (engine stopped), L 2
if x 2 = 0 and u 3 = 1 (engine start-up).

(3)

While x 1 (t) is a continuous state governed by ẋ1 (t) = f 1 (u(t), t), the state x 2 (t) is a discrete state, and takes only the values 0 or 1 corresponding to engine stopped/engine started. The introduction of x 2 (t) permits to take into account an extra fuel consumption, L 2 , corresponding to the necessary energy to start the engine. Physically, this energy comes from the inertia of the engine, and is consumed by the starter (electric energy), and by the engine (small quantity of fuel to reach the idle speed). The value L 2 is considered as a constant.

The general optimization problem is the following:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ min u∈U J (u) := T 0 L(x(t), u(t), t)dt + g(x 1 (T ), T )
subject to:

ẋ1 (t) = f 1 (u 1 (t), t), x 1 (0) = x 1,0 , x 2 (t + ) = x 2 (t) + u 3 (t), x 2 (0) = x 2,0 x min 1 ≤ x 1 (t) ≤ x max 1 (4)
where 0 and T are respectively the initial and the final times of the prescribed driving cycle, and f 1 represents the dynamics of the battery. This minimization problem must be solved on a prescribed cycle, for instance the NEDC cycle, Fig. 1. Thanks to a vehicle model (mass, wheel diameter, aerodynamic coefficients, etc.) and the vehicle speed request V v (t), the requested torque T w (t) and speed ω w (t) at the wheels are computed. In this optimization problem, we consider the cycle as a wheel torque T w (t) and wheel speed ω w (t) trajectory to be followed.

System modelling

In this application, we consider that the engine is located upstream of the gearbox, as on conventional vehicles. Besides, the electric motor is located downstream of the gearbox, beside the wheels. Let us introduce now some useful notations to describe the powertrain and to define the control u(t):

• the requested torque at wheels is written as T w (t) = R f ar R gr T e (t) + T m (t) + T bk (t), and can be achieved thanks to the engine torque T e (t) multiplied by the front axle ratio R f ar and by the gear ratio R gr , and thanks to the motor torque T m (t). A brake torque T bk (t) can also be applied by braking system.

• the engine torque is T e (t) = u 1 (t)T w (t)/(R f ar R gr ), • the electric motor torque is T m (t) = (1 -u 1 (t))T w (t).
During braking phases, T rq (t) < 0, the braking torque can complete the electric motor torque to slow down the vehicle: T bk (t) = T w (t) -T min m (t) with T min m (t) the minimal electric motor torque depending on the engine speed, thus depending on time.

A simple model of the battery is implemented, the state of charge of the battery evolving as

ẋ1 (t) = f (u 1 (t), t) = - ω m (t)T m (t)K U batt (t)n capa (5)
with ω m (t) the electric motor speed, U batt (t) the battery voltage, K an efficiency, and n capa the nominal capacity of the battery. In the following, we assume that the battery voltage is constant along the cycle. We set K = K U batt .n capa and by replacing T m (t) by its expression, we obtain

ẋ1 (t) = -Kω m (t)(1 -u 1 (t))T w (t). ( 6
)
This very simple model for battery dynamics may be complexified to take into account resistive loses.

The state equation for x 2 is ẋ2 =0. It is easier to consider x 2 in discrete time to understand its behavior (see Sect. 3).

In this optimal control problem, we make several assumptions about the powertrain: (i) the wheel speed ω w is equal to the electric motor speed ω m , (ii) as the optimization problem is applied on a full hybrid vehicle, the engine can be stopped when the torque is provided only by the electric motor. As mentioned above in the definition of L(x(t), u(t), t), a fuel consumption penalty corresponding to the necessary energy to start the engine must be considered, (iii) recharging the battery is possible for negative torques (braking request) or by an additional engine torque beyond the driver request torque. In the last case, the command u 1 (t) is larger than 1, (iv) we consider that the average efficiency of the electric motor is 0.85.

As in the general optimal control problem (4), the command variables u(t) and the state variables x(t) are submitted to several constraints:

• the engine can only produce a positive torque, and is limited to a maximal torque which depends on engine speed ω(t): 0 ≤ T e (t) ≤ T max e (ω(t)), • the electric motor torque is bounded:

T min m (ω(t)) ≤ T m (t) ≤ T max m (ω(t)),
• the storage capacity implies a minimal and a maximal state of charge of the battery (which are fixed to 50% and 70% in our study):

x (1) min ≤ x 1 (t) ≤ x (1) max ,
• the final state of charge is constrained to be close to the initial state of charge to maintain a null electrical energy balance (to avoid to discharge totally the battery for minimizing the consumption) x 1 (0) = x 1 (T ).

Notice that T w (t) and ω w (t), which depend on the vehicle model, are computed thanks to a quasistatic approach [START_REF] Beck | Comparison of two real-time predictive strategies for the optimal energy management of a hybrid electric vehicle[END_REF][START_REF] Tate | Finding ultimate limits of performance for hybrid electric vehicles[END_REF][START_REF] Anatone | Control-oriented modeling and fuel optimal control of a series hybrid bus[END_REF]:

ω w (t) = V v (t) r wh , (7) T w (t) = - 1 r wh mg sin α + Fr(V v (t)) + 1 2 ρSC x (V v (t) + v 0 ) 2 -m + J wh r 2 wh dV v (t) dt , ( 8
)
where r wh is the wheel radius, α is the road slope, 1 2 ρSC x (V v (t) + v 0 ) 2 corresponds to the aerodynamic loses, Fr(V v (t)) is the rolling resistance force, m is the vehicle mass, J wh is the wheel inertia, and v 0 corresponds to the wind speed (null here).

Constraints on command u 1 (t) are directly derived from the constraints on maximal and minimal engine and electric motor torque

u 1 (t)T rq (t) ≤ T max e (t), (9)
T min m (t) ≤ (1 -u 1 (t))T rq (t) ≤ T max m (t). ( 10
)

Dynamic programming

The Dynamic Programming method (DP) is classically applied to solve the problem (4) [START_REF] Scordia | Systematic elaboration of online energy management laws for hybrid vehicles[END_REF][START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF]: it relies on a very simple idea, the principle of optimality, stated by Richard Bellman. It states that an optimal policy can be built step by step, first by building an optimal control for the "tail subproblem" involving the last step, then by extending the optimal policy to the "tail subproblem" involving the last two stages, and by iterating until an optimal control has been built for the complete time interval. Thus, the optimization problem is solved backward from final time to initial time. The criterion to be minimized in (4) is thus discretized in time, leading to

min u∈U J (u) := N -1 k=0 L(x k , u k , k) t + g(x 1,N ), ( 11
)
where L(x k , u k , k) t is the cumulated fuel consumption over the time interval

[k, k + 1],
x k is the state vector at step time k, and N is the final time of the driving cycle. The constraint on the final state of charge x 1,0 = x 1,N has been introduced via a penalization term g(x 1,N ) = β(x 1,Nx 1,0 ) 2 , with β, a constant to be chosen.

In general, β is most probably different for battery size. However, the same value of β, chosen high enough, has been used for all the tested battery capacities.

Notice that this minimization must be done under state constraints x 1 ≤ x max and x 2 ≥ x min , and the control

(u 1 , u 2 , u 3 ) = u ∈ U = {U 1 , U 2 , U 3 } with ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ U 1 (x 2 , k) = [u min 1 (x 2 , k), u max 1 (x 2 , k)], U 2 (k) = tg min (k), . . . , g max (k)u, U 3 (x 2 ) = {-1, 0} if x 2 = 1, {0, 1} if x 2 = 0, (12)
where g min (k) and g max (k) correspond to minimal and maximal allowed gears (which are discrete values) at time step k, and u min 1 (x 2 , k) and u max 1 (x 2 , k) come from ( 9) and (10) expressed at time step k.

From Bellman principle, the minimum cost-to-go V k (x k ) for the time interval

[k, N], 0 ≤ k ≤ N -1, is then expressed as follows V k (x k ) = min u k ∈U L(x k , u k , k) t + V k+1 (x k + f (x k , u k , k) t) , ( 13
)
where f (x k , u k , k) is the system dynamics given by:

x k + f (x k , u k , k) = x k+1 = x 1,k+1 = x 1,k + f 1 (u k , k) t, x 2,k+1 = x 2,k + u 3,k . ( 14
)
At time N , the cost function is

V N (x N ) = g(x 1,N ). ( 15
)
This optimization problem is solved backward from final time step to initial time step using a discretization of function V in the control space and in the state space, that gives a 3-D grid (time t, state of charge x 1 , and engine state x 2 nodes). A standard time step used in our examples is t = 0.5 s, the discretization step x 1 for the state x 1 is 0.25%, and x 2 takes only 0 or 1 values. Thus, the number of state nodes is M = 2((x max 1 x min 1 )/ x 1 ). A standard DP algorithm consists in determining, at each step k and at each state value, the control feasible set U(x k , k), and calculating the optimal control that minimizes the sum

(L(x k , u k , k) t + V k+1 (x k + f (x k , u k , k) t)).
Once initial time step is reached, a control law synthesis can be done, starting from initial conditions, minimizing again L(x k , u k , k) t +V k+1 (x k +f (x k , u k , k) t), and integrating state equations.

As the number of control is relatively small, about one hundred chosen controls, it's easy to select the optimal one by a simple minimization procedure.

This method allows a fast optimization: (N -1)M constrained optimization problems (13) of only three parameters u 1,k , u 2,k and u 3,k are solved with a fine time discretization of the variables. This provides the optimal power repartition u * k at each time step (we refer to [START_REF] Guilbaud | Control optimization for a power-split hybrid vehicle[END_REF] for some theoretical results on the convergence of the method and error estimations).

A direct constrained optimization of ( 11) is an alternative method but leads to a large non linear optimization problem with a large number of inequality constraints. This is the reason why the DP algorithm is often preferred when the number of state variables is small (one or two state variables). 4 Variations of the size of some powertrain elements

On conventional vehicles, the maximal engine power is mainly defined by its displacement. Economic vehicles have a small engine, with a small fuel consumption, but with a low maximal torque. When engine displacement is higher, fuel consumption is higher, but maximal engine torque is also higher.

In hybrid vehicles, engine displacement has still almost the same impact on fuel consumption. However, the presence of a second energy source, which implies at least an electric motor and a battery, gives other degrees of freedom to decrease the fuel consumption.

In this study, the optimization of gear shifting is not presented, because gear shifting is imposed on the NEDC cycle. However, it is possible to impose a penalty at each time the gear is shifted.

The following results have been obtained by testing different sizes of elements of the powertrain. This approach has been preferred to parametric optimization, which is a quite complicated problem. The reader can nevertheless see some examples in [START_REF] Bonnans | Parametric optimization of hybrid car engines[END_REF].

Reference case

The reference case is a Citroën C1 with a 1000 cc gasoline engine. We consider on this vehicle that an electric motor of 10 kW is located downstream of the gearbox, at the same level as the wheels, and that a battery of 1.3 kW/h provides the electric energy to the motor. The potential of this vehicle is assessed on the NEDC cycle, and gives the results summarized in Table 1.

Variations of battery capacity with constant battery weight

The battery dynamic (5) depends on the battery capacity: for a constant control, the battery dynamic is bigger as the battery capacity is smaller. As a result, the state of charge limits, x min 1 and x max 1 , can be more often reached with a very small battery than with a big battery.

As the system dynamic does not depend on the state x 1 , the fuel consumption will not evolve with the battery capacity as soon as the state of charge remains inside its bounds x min 1 , x max 1 without touching them. Table 2 gives optimal fuel consumption on NEDC cycle for varying capacities.

Variations of battery capacity and battery weight

By considering the variations of battery capacity and the corresponding variations of its weight, not only the available electric power is changed, but also the global weight of the vehicle. In this context, the necessary power to match the requested vehicle speed of the cycle becomes higher as the weight of the battery increases. Thus, the variation of battery capacity causes the engine to change its operating points: same engine speeds, but different requested torques. Also, the vehicle will potentially be able to recover much more electric energy during braking, if control constraints are not activated. Consequently, the analysis of these results is not trivial. We consider in the following that the battery weight varies according to 1 kg for 50 W h. The Table 3 shows obtained results.

It is noticeable that those results show a different trend than the results of Table 2. An optimal value is reached for n capa = 1.04 kW h, as seen on Fig. 2.

Variations of maximal and minimal electric motor torques

The car manufacturers are interested in estimating the best size of an electric motor for a hybrid vehicle, according to its use, as there is no need to have a very powerful motor if half of its maximum power is never reached.

The variations of maximal and minimal electric motor power may change the fuel consumption, but it also changes the state of charge trajectory (see results in Table 4).

For this vehicle used on the NEDC cycle, the fuel consumption does not evolve anymore if the electric motor power is larger than P = 15 kW. Smaller powers of electric motor do not affect much the fuel consumption, as the gain for an electric motor is still 20.15% for a power of 5 kW. 5 Real-time control strategy

Validation of DP results on a more realistic model

Results obtained with the dynamic programming applied on the Matlab model are tested on a more realistic model implemented in AMESim, corresponding to the reference case, see Fig. 3. By comparison with the Matlab model used for optimization, this AMESim model includes all axle inertias, lag times of engine and electric motor, and a driver model which controls acceleration and brake pedal to follow the speed request from the cycle. Figure 4 shows the behavior of AMESim model with the optimal control of electric motor obtained from dynamic programming on the Matlab model. Figure 4(a) shows the effective vehicle speed (blue), compared with the speed request from the cycle (red). Figure 4(b) and(d) show the motor and engine torques, while Fig. 4(e) specifies if the engine is started or stopped. Figure 4(c) shows a comparison between optimal SOC trajectory obtained with Matlab model, and the one obtained with AMESim model.

The global fuel consumption obtained with AMESim model is 310.2 g (3.81 l/100 km) while the one with Matlab model is 310.34 g (3.75 l/100 km), but the vehicle linear displacement is not the same: Matlab model covers 11.01 km, but AMESim model only covers 10.84 km, because of inertia that cause a bad follow-up of the vehicle cycle at low speeds. The resulting error between the two models is about 1.5%. In this section, suboptimal control law is derived from the optimization results. The chosen real-time control law is the Equivalent Consumption Minimization Strategy [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF]Guzzella and Sciarretta 2005;[START_REF] Paganelli | Equivalent consumption minimization strategy[END_REF]: it is based on instantaneous equivalent fuel consumption, the battery being considered as a second energy source. ECMS aims to choose the cheapest energy source, between fuel tank and battery storage.

This strategy is inspired by the Pontryagin's Principle, which set optimality conditions for the unconstrained continuous optimization problem (4). Some important theoretical results can be found in [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]. To apply this strategy, we form the Hamiltonian function H, associated to the problem (4) without considering the state constraints, and given by

H(x(t), u(t), p(t), t) = L(x(t), u(t), t) + p(t) ẋ1 (t), ( 16
)
where L(x(t), u(t), t) t represents the variation of fuel consumption, and given by

L(x(t), u(t), t) = L 1 (u(t), t) if x 2 = 0, 0 i fx 2 = 0, ( 17
)
ẋ1 corresponds to the variation of state of charge over a time step, and p is a Lagrange multiplier associated with the constraint ẋ1 (t) = f 1 (u(t), t). In this strategy, the optimal control is the admissible one which minimizes H, taking into account all the control constraints. The same models as the ones used in dynamic programming algorithm are used in ECMS, which is the advantage of this model-based strategy.

Fig. 4 Results of optimal control applied on AMESim model

Meaning of p

The coefficient p corresponds to the price of electric energy. The larger is p, the more expensive is the electrical energy, and the more interesting to be recovered. The smaller is p, the cheaper is the electrical energy, and the more interesting to be used to power the vehicle conjointly with the engine. The variation of p depends on the time evolution of the dynamic state variable x 1 , according to the Euler-Lagrange equations, namely

ṗ(t) = - ∂H ∂x 1 (x(t), u * (t), p(t), t). ( 18
)
As either L or ẋ1 do not depend on x 1 , (18) gives

ṗ(t) = 0. ( 19
)
As a results, when no state constraint is active, p is constant for the considered system (6). When a state constraint becomes active, the sign of ṗ is completely known, depending on the active constraint: x = x min or x = x max . Indeed, if the augmented Hamiltonian is written, with state constraints, we get

H(x(t), u(t), p(t), t) = L(x(t), u(t), t) + p(t) ẋ1 (t) + λ 1 (t)(x min 1 -x 1 (t)) + λ 2 (t)(x 1 (t) -x max 1 ), (20)
where λ 1 ≥ 0 and λ 2 ≥ 0 are Lagrange multipliers associated with state constraints x min 1x 1 ≤ 0 and x 1x max 1 ≤ 0. Solving ṗ = -∂H ∂x 1 leads to ṗ = λ 2λ 1 . One can verify that if x 1 = x min 1 , then λ 2 = 0 and λ 1 ≥ 0 meaning that ṗ ≥ 0. The same idea is used with the second state constraint.

Control of p

In ECMS, the value of p is often controlled thanks to rules of [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF], Liu and Peng (2006), mainly depending of the state of charge. Indeed, as the future is not known, it can be better to keep the state of charge x 1 far from its bounds, and to control p as a function of x 1 with a feedback controller. However, in this paper, we consider that p remains constant, and we just try to show the potential of this strategy with an adequate initial value of p.

An estimation of p can be obtained by calculating the equivalent Lagrange multiplier with the dynamic programming results. From [START_REF] Bryson | Applied optimal control[END_REF], equivalence between Pontryagin's principle and dynamic programming in continuous time can be found, assuming that V is differentiable, with

p(t) = ∂V ∂x 1 (t, x 1 ), (21)
where V (t, x 1 ) corresponds to the cost-to-go function.

To estimate the value of p, we calculate the partial derivative of V (defined in ( 13)) with respect to x 1 with finite differences, for each step k. Using theses results allows to estimate that the average value1 of p should be into the range [-3400; -3000]. Because of some differences between the model used with dynamic programming and the AMESim model, a trial & error method has been applied to end up with p = -3075.2 This value allows to reach a final state x 1 (T ) close to the initial state x 1 (0).

Results of the ECMS Real-time controller

The Fig. 5 shows some results obtained with this real-time control strategy, with a correct value of p.

Figure 5(a) shows that the hybrid vehicle still follows the speed request from the cycle. However, some differences appear on the state of charge, Fig. 5(b). The fuel Fig. 5 Results obtained with a sub-optimal strategy applied on AMESim model consumption obtained with this sub-optimal strategy is 326.8 g, corresponding to 3.97 l/100 km. Compared to conventional vehicle, this controller-with an ideal value of p in this example-allows to save 16.6% of fuel on the NEDC cycle. Some rules to control p with the state of charge should be implemented to make this controller more robust.

Conclusion

In this paper, has been presented a study of a full hybrid vehicle, with a classical offline optimization tool, and an example of a real-time control strategy inspired from optimal control theory.

A simple model in Matlab has been implemented and coupled with a dynamic programming algorithm. Through this model, an engine startup cost L 2 has been introduced in the optimal control problem, thanks to a discrete variable x 2 , to take into account the extra fuel consumption of the engine in order to reach its idle speed. This penalization term in the cost function avoids too many successive engine stops and engine starts, and is usually never considered in this kind of optimal control problems.

This optimization tool allows to estimate the optimal sizes and powers of the elements of hybrid powertrain. Nevertheless, these results, obtained on the NEDC cycle, depend on the chosen driving cycle: so the optimal sizes can be a little different on other cycles.

After a validation of the optimal trajectory and torques applied on a more realistic model in AMESim, a model-based sub-optimal controller has been implemented, inspired by Pontryagin's Principle and optimality conditions. In the paper, we briefly show how to take into account some dynamic programming results to initialize the value of p.

However, even if an average value of p can be estimated from offline optimization results, which provides an initial value, the choice of an adequate initial value remains a delicate task, and p should be tuned regarding to the evolution of the state of charge of the battery, and to the type of the cycle.

Fig. 1

 1 Fig. 1 NEDC driving cycle: vehicle speed at left, gear ratio at right

Fig. 2

 2 Fig. 2 Variation of gain of fuel consumption with battery capacity

Fig. 3

 3 Fig. 3 Validation AMESim model

Table 1

 1

	Fuel consumption for the reference case	Fuel consumption (l/100 km)	Gain (%)
	Conventional vehicle	4.76
	HEV	3.76	20.99%

Table 2

 2 Study of the impact of

	the battery capacity		(l/100 km) Fuel consumption	Gain (%)
	Conventional vehicle		4.76
	HEV with battery	0.325 kWh 3.78	20.40%
		0.65 kWh	3.76	20.99%
		1.3 kWh	3.76	21%
		2.6 kWh	3.76	21%
		5.2 kWh	3.76	21%
	Table 3 Study of the impact of variations the battery capacity with weight		(l/100 km) Fuel consumption	Gain (%)
	Conventional vehicle		4.76
	HEV with battery	0.325 kWh 3.78	20.40%
		0.65 kWh	3.76	20.84%
		1.04 kWh	3.76	21.01%
		1.3 kWh	3.76	20.99%
		1.625 kWh 3.76	20.90%
		1.95 kWh	3.76	20.79%

Table 4

 4 Study of the impact of the maximal motor power

		Fuel consumption	Gain (%)
		(l/100 km)
	Conventional vehicle	4.76
	HEV with motor P = 5 kW	3.80	20.15%
	HEV with motor P = 10 kW	3.76	20.99%
	HEV with motor P = 15 kW	3.74	21.25%
	HEV with motor P = 20 kW	3.74	21.25%

In practice, a shooting algorithm can be used by setting the initial value of p, integrating the system equations, and iterating on the initial value of p to reach the desired final state of charge.

If the variations of fuel and chemical energies are used instead of the variations L(u(t), t) t and ẋ1 (t) t , the coefficient p can be scaled to a smaller value.


D Sinoquet

Modeling a Priori Information on the Velocity Field in Reflection Tomography

Reflection tomography consists in determining a velocity field and reflector geometries from traveltimes picked on multioffset seismic section. The solution of the tomographic inverse problem being underdetermined, we need to integrate a priori information on the model. A regularization by means of model curvature has in general no physical justification and leads to geologically incorrect models. This paper presents a formulation that integrates in a realistic way a priori geological information associated with the regularity of the model and with the relation between the velocity distribution and the interface geometries. Such an integration seems particularly critical when dealing with smooth velocity models in which velocity and interfaces are defined independently from each other. We validate the interest of our formulation on a real data example.

INTRODUCTION

Conventional seismic processing can not image complex geological structures such as the one shown in Figure 1. Lateral velocity variations in the vicinity of the fault prevent complete imaging of the base of salt and of the underlying structure.

Successful seismic imaging of such geologic structures requires prestack depth migration which itself requires the determination of an accurate velocity model.

Reflection tomography (Bishop et al., 1985, Chiu et al., 1987) seems very attractive for this determination. In this method the velocity distribution and the geometry of the interfaces are determined from traveltimes picked on prestack seismic sections.

Because of the finite number of data and of the infinite number of unknowns the tomographic inverse problem is ill posed. To offset the indetermination additional information on the model is necessary. Delprat-Jannaud and Lailly (1993) proposed a regularization by means of the model curvature that insures the uniqueness of the solution and mathematical stability. Unfortunately this method is unstable in practice when applied to complex geologic structures where the curvature of the velocity can, in no way, be understood as smooth. As a consequence of this unstability the velocity distribution and the structural model can be found completly uncorrelated from each other as shown in Figure 2. Even if the model matches the picked traveltimes, it is geologically unacceptable.

REFLECTION TOMOGRAPHY ON COMPLEX GEOLOGIC STRUCTURES

SMOOTH VELOCITY MODELS

We are primarly interested in applications of reflection tomography for the determination of a velocity model to be used as input to prestack depth migration to image complex structures. Discontinuous velocity models are not suitable for such a migration: misplaced discontinuities in the velocity distribution create discontinuous slopes in the migrated reflections that can completely mislead the interpreter. To preserve the nice features of reflection tomography, namely its ability to produce a velocity model kinematically consistent with the seismic data, the velocity model used for the migration should be identical to the one obtained by tomography. This lead us to build a reflection tomography dealing with smooth velocity models (Figure 2).

TOMOGRAPHY PROBLEM

Reflection tomography involves the solution of two problems: the forward problem and the inverse problem. Forward problem: We deal with two-dimensional media limited to a rectangular domain defined by a horizontal interval and by a depth interval The model m is composed of a slowness squared distribution that is supposed to be and of explicit interface depth functions supposed to be i indicating the number of interface ( i = 1, . . . , N ). The forward problem is the computation of traveltimes in a given model by raytracing. We consider acoustic waves propagating away from a source, reflecting on an interface (reflector) and propagating to a receiver at the surface. A ray R j is defined by the triple (source, reflector, receiver) indexed by j. We define the corresponding traveltime as Inverse problem: We have at our disposal a set of observed traveltimes T obs and we look for the model that matches these data at best. The least square formulation of the inverse problem consists in minimizing the misfit between observed and computed traveltimes : More generally, the solution of reflection tomography is basically underdetermined and relevant geologic information has to be integrated to remove this underdetermination. This underdetermination is quite critical in complex structures where difficulties in picking seismic reflections make the traveltime data very sparse. In this paper we propose a method for the integration of some specific geologic information which is particularly suited for applications to complex structures.

where is a norm in data space, for instance the Euclidian n o r m . This non linear inverse problem is solved by the Gauss-Newton method which consists in solving successive linearized problems that can be expressed as the minimization of the functional ( 2) with a given model, the Jacobian of the forward modelling operator , a model perturbation giving the new model These linearized inverse problems are underdetermined: they admit several solutions. Additional information is necessary to overcome this difficulty.

A PRIORI INFORMATION ON VELOCITY STRUCTURE

The subsurfaces of interest are not chaotic media: the interpreter knows the sedimentary deposition modes and the unconformity surfaces from regional geology. A certain velocity structure results from such an organization. We use two types of geological information on the velocity field: information on the regularity of the velocity and on the anisotropy of regularity due to sedimentation. (Sinoquet and Brac, 1992) Regularity of the velocity: For complex structure tomography, we deal only with macroscopic phenomena: very slight velocity variations will not be determined by traveltime inversion. Versteeg (1991) showed that only velocity field wavelengths greater than 100 meters are of importance for prestack migration. Accordingly, this defines the objective for complex structure tomography. Additionally, geology tells us that the velocity is more regular in some places than in others, even though it is assumed smooth. Figure 3 shows the wide velocity variations in the medium, particularly in the neighborhood of the fault and at the base of the salt, as well as at the hard and soft carbonate contacts (at the right of the fault) and between these carbonates and the clay-sand sequence. Anisotropy of regularity (influence of sedimentation):

The sedimentation makes regular, roughly horizontal deposits. In turn the velocity varies much less along deposition lines in the orthogonal direction (Figure 4). Therefore the velocity varies little along the reflectors that are geological isochrons.

To model the above described a priori information we add penalization terms to the tomography objective function

(3) where is the gradient of the slowness squared and is the unit tangent vector to the interface The velocity guide term constrains the slowness squared (consequently the velocity) to vary little along certain reflectors: it integrates the second type of a priori information. The regularization terms constrain the model to be regular and help to propagate the local velocity guide on the interfaces into the layer. But the velocity regularization term plays also an active role in making the velocity distribution in accordance with the reflector geometries.

By imposing gentle variations of the velocity gradient, we make this gradient remain approximately orthogonal to the geological isochrons that are within a layer between the reflectors involved in the application of reflection tomography. Consequently the isovelocity lines and the geological isochrons are expected to have similar directions within a whole layer.

The variable weighting functions represent the uncertainty associated with each piece of information. They enable us to calibrate the functional terms according to their physical meaning. They are space varying and depend also on the model to allow for instance possible degenerescence of the velocity regularization in the vicinity of faults whose geometry is an unknown of reflection tomography. We minimize the non quadratic objective function (3) by the Gauss-Newton method. Each linearized problem is solved by a block relaxation method: we update velocity parameters and interface parameters alternatively and we iterate the process. Such a method avoids ill conditionning difficulties involved in the minimization of a function that depends on parameters of different physical natures. At each relaxation step we use an algorithm of constrained optimization based on the Augmented Lagrangian method, developped by Glowenski and Tran (1993). It allows to integrate linear constraints such as imposing an interface slope or imposing the velocity to be smaller than a certain value.

COMPARISON OF MODELS OBTAINED BY TOMOGRAPHY WITH DIFFERENT A PRIORI INFORMATION

The model in Figure 2 was obtained by tomography with curvature regularization (Delprat-Jannaud and Lailly, 1993), whereas the model in Figure 5 is the result of tomography integrating the a priori information described above. Both models match the traveltime data.

Basically, models associated with complex structures do not have small curvature. Thus to obtain a model that matches the traveltime data we can only use very weak regularization, thus leading to unstability in pratice (Delprat-Jannaud and Lailly, 1993). Therefore it is not surprising that the model in Figure 2 is geologically unacceptable. Indeed, there is no coherency between the velocity field and the interface geometries. Moreover, as a consequence of the weak curvature regularization, unrealistic velocity anomalies appear: they reflect the underlying unstability.

To obtain the model in Figure 5, we have constrained the velocity to vary little along the interfaces, except along the fault and the base of salt. All the interfaces has been regularized by second derivatives. Also the second derivatives of the velocity are penalized everywhere but with different weights according to the region. Indeed, we define three zones: the neighborhood of the fault, where the regularization weight is very low to allow wide velocity variations, and the two regions on both sides of the fault.

All the unrealistic features of the model in Figure 2 have disappeared in the model in Figure 5: • the velocity varies regularly within each layer • the isovelocity lines follow nicely the interfaces.

In the neighborhood of the fault, the velocity is not well determined because of the weak regularization introduced in this area.

So, we obtained a model which matches the traveltimes and which is geologically acceptable. But as we will see in the next section, a difficulty remains to obtain such a model: the choice of the weights associated with the a priori information.

INCIDENCE OF THE CHOICE OF WEIGHTS

The model in Figure 6 was obtained in the same way as the model in Figure 5 but with different weights: we increased weights associated with the velocity regularization and with the velocity guide.

This model is very regular on both sides of the fault and presents high velocity variations at the top and the bottom of the fault. The model respects the a priori information but does not match traveltimes. This means that the information we try to integrate is not consistent. Therefore a good compromise has to be found between the conflictual pieces of information and the choice of weights has to be adequate.

In fact, we have to choose their values according to the physical uncertainties they represent. This is far from obvious a priori. However, an iterative trial and error approach turned out to be quite constructive to calibrate those weights. If the weights guessed a priori have led to the model in Figure 6 we conclude that, the traveltime match being unacceptable, we should increase the associated weight or, equivalently, decrease the weights associated with the velocity guide or with the velocity curvature or both. This trial and error approach was effectively used to calibrate iteratively the weights so as to obtain the model in Figure 5.

C O N C L U S I O N

Integration of a priori geological information especially on the smoothness of the subsurface is crucial in reflection tomography: it is the warrant of stability.

Integrating information ensuring coherency between the velocity distribution and geometry of the interfaces is essential when tomography makes use of a velocity distribution represented by a continuous function. Such a representation is appropriate when the velocity computed by reflection tomography is to be used as input to prestack depth migration for the imaging of complex structures. Thus we proposed a modeling of the regional regularity of the velocity and of the sedimentation effects on its structure.

Solving a tomography problem integrating such a priori information leads to the minimization of a multi-objective function. A weight is associated with each term of the functional. These weights have been model space varying to allow degenerescence of the velocity regularization in the vicinity of faults, for example. The tuning of these weights is quite critical because it defines the compromise to choose between conflictual pieces of information.

A correct tuning has been obtained by a trial and error approach, thus leading to a model matching both traveltime data and the a priori geologic information.

  The linearized problem associated with the minimization of the objective function is a well posed problem: as in Delprat-Jannaud and Lailly (1993) uniqueness and stability of the solution are ensured. NUMERICAL RESULTS NUMERICAL METHOD Parameterization of the model: The model is represented by cubic B spline functions, adequate to build functions (Barsky et al., 1987). The velocity parameters are the vertices of a tensor of B spline functions (respectively depending on and on The interface depths are represented by B spline function depending on Optimization:

ACKNOWLEDGMENTS

This research was carried out as part of the Prestack Structural Interpretation consortium project (PSI). I hereby acknowledge the support provided by the sponsors of this project.

I also wish to thank Patrick Lailly, Jean Brac and Florence Delprat-Jannaud for their valuable advice. I want to thank Elf Aquitaine for its financial support.

R E F E R E N C E S

Barsky, R.H., Bartels, B.A. and Beatty, J.C., 1986, An introduction to splines for use in computer graphics and geometry modeling, Morgan Kaufmann Publishers Bishop, T.N., Bube, K.P., Cutler, R.T., Love, P.L., Resnick, J.R., Shuey, R.T., Spindler, D.A., and Wyld, H. W., 1985, Tomographic determination of velocity and depth in laterally varying media, Geophysics, vol.50, 903-923 Chiu, S.K.L., and Stewart, R.R., 1987, Tomographic determination of three-dimensionnal seismic velocity structure using well logs, vertical seismic profiles, and surface seismic data, Geophysics, vol. 52, 10854098 Delprat-Jannaud, F. and Lailly, P., 1993, Ill posed and well posed formulations of the reflection tomography problem, accepted for publication in J.G.R. Sinoquet, D. and Brac, J., 1992, Modeling a priori information on velocity field for reflection tomography, PSI 1992 Annual Report, Institut Francais du Petrole, Rueil, France Glowinski, R. and Tran, Q.H., 1993, Constrained optimization in reflection tomography: the Augmented Lagrangian method, submitted to Journal of optimization, SIAM Versteeg, R.J., Analyse du probleme de la determination du modele de vitesse pour l'imagerie sismique, these de Doctorat, Universite Paris VII, France


domain: info sdu

text:
SUMMARY

The SMART (Sequential Migration Aided Reflection Tomography) method, as explained in the first part of this paper, starts after a first set of traveltimes in the unmigrated prestack data has been picked and the inventarization of useful a priori knowledge related to these traveltimes has been made. Thereto a preparative phase is needed. First a global estimate of the subsurface structure is made. Hereto we use the standard stacking and poststack interpretation procedures which 'allow for getting insight in the degree of complexity of the subsurface. Next the traveltimes can be picked. When interpreting prestack data important qualitative structural information in difficult target zones (e-g. fault zones or salt structure flanks) can be obtained. Such an analysis guides the interpreter in selecting and picking the best traveltimes of primary events. Once the preparation is finished the SMART method can be applied for a detailed determination of a structural and velocity model in a very consistent way. It is emphasized that velocity variations in complex structures can be determined accurately by prestack traveltime inversion techniques. This phase has an iterative character. In order to update the velocity model after the first iteration additional traveltimes are needed. Next additional traveltimes are obtained by interpretation of the cube of migrated data which can be easier than in the time domain due to the focussing and positioning effect of the migration process. By tracing rays in the same velocity model as was used for mi.gration on the newly interpreted events, we will obtain additional traveltimes which will make the set of input data for the next iteration of tomography more complete. A new velocity model is calculated and the data are remigrated. In this paper we will demonstrate the feasibility of this approach using a 2D real data set. We executed a number of iterations of the SMART method and ended up with of the complex structure. a very satisfactory depth image

THE DATA

We used for this application a 2D dataset covering a salt structure. It consists of 300 shotrecords at a regular interval of 40m. The acquisition was done in a split spread. The half spread length is 1920 meters with 48 geophones. The data were delivered with a standard preprocessing (filtering, zero-phase deconvolution and muting). Because of some clearly visible groundroll, we applied a second filter in order to remove most of this in Figure 1. low frequency noise. A partial stack of the data is shown

THE PREPARATIVE PHASE

Analysis of complexity

In order to get an idea of the degree of complexity of a subsurface, it is useful to construct several partial stacks with the same stacking velocity model. Because the stacking process is based on flattening of the hyperbola's in CMP's, through some NMO and DMO based correction, differences in between the partial stacks demonstrate the failure of the process. In areas with complex subsurface structures these hyperbola's aren't necessarily flat due to different raypaths left and right of the midpoint. In this dataset this phenomenon can be observed in a series of CMP's covering the saltdome (See Figure 2). Another way to get an idea of the complexity is to do a post stack depth migration by a layer stripping approach using the best partial stack. For these data the results are satisfactory for the sedimentary zones left and right of the dome, but are incorrect for the deep interfaces and the base of the salt. This is partially due to events that are lost during the stacking procedure. Other causes for this failure are: the uncertainty in picking the right interface that serves as the next velocity boundary and the difficult choice of the velocities which becomes more and more hazardous as the depth increases. The final result is unreliable and the resulting depth for the base of the salt depends largely on the choices made by the interpreter Clearly these data cannot be handled by standard processing techniques.

Left and right of the salt dome and below it the nature of the trace gathers is too complex. A prestack imaging method using a velocity model computed by tomography seems adequate for solving the aforementioned problems.

Data preparation for the SMART method

The next step after the analysis of the complexity is the data preparation for the SMART method. Its goal is to prepare an initial set of traveltimes to be used in the first iteration. We split this phase in a number of consecutive sub-phases:

• Creating a initial set of guides for the prestack interpretation.

• Picking traveltimes.

• Quality control of the traveltimes.

• Selection of representative traveltimes and calculation of the associated weights.

Creating a set of guides.

Guides are indicators for the interpreter suggesting where to look in the prestack unmigrated data for a certain event. They are also warnings for complicated situations as multiples, triplications and situations were no reliable indications for the nature of an event is available. The geologic guides are qualitative (e.g. presence of a fault) or quantitative (e.g. the depth of horizon A is 2500m). The geophysical guides are for example the presence of multiples or diffractions. They are derived from the unstacked or stacked data.

For this dataset the following data were used: a set of (partial) stacks, time-and depth-migrated stacks and the cube of preprocessed prestack data. It allowed us to determine the zones where picking traveltimes directly in the unmigrated data could lead to incorrect traveltime information for the tomography. These zones are indicated in Figure 1 (Za and Zb, a zone with triplications and a series of unexplained events.

Picking the first set of traveltimes

Using the guides the picking of the traveltimes can start. This is done in the cube of unmigrated data. There is no preference for picking in a specific trace gather. This depends of the available guide. When it is a geological one the common offset gathers are most suited. Using a geophysical one the interpretation is done in the shotgathers or the common midpoint gathers. Whatever direction is chosen, one has to end up with one consistent surface describing the set of prestack traveltimes for one reflector.

Quality control of the traveltimes.

For those reflectors for which split spread data are available we analyzed to what extend the principle of reciprocity is fulfilled. The criterion of reciprocity can help the interpreter to find gross errors in his interpretation. An example is shown in Figure 3. There is a strong correlation visible between traveltimes along one specific offset. This is an indication that the offset has been picked incorrectly. The interpretation has been corrected and the strong offset-related errors were removed.

A second control consists of the calculation of the difference between an ideal hyperbola and the picked traveltimes. Given the traveltimes sorted per CMP, the ' stacking velocity' for each CMP can be calculated. It informs the interpreter about the regularity of his interpretation. We call this the CMP difference plot. Going from one geophone position to another the difference in traveltimes should be reasonable; that is to say within the limits of the standard interpretation error which is usually estimated to be 8 ms.

Selection of representative traveltimes and calculation of the associated weights. Not all of the traveltimes picked are really needed for our implementation of reflection tomography. It is possible to give weights to each picked traveltime so that in less complex parts of the seismic data a selection of traveltimes is sufficient for obtaining the same result as if all traveltimes were used, of course at a significant lower computational cost. The selection of representative traveltimes is based on the analysis of the aforementioned CMP difference plots. They give a qualitative measure of the regularity of the velocity field above the reflector. For these data it was decided in the less complex parts to resample the set of traveltimes at every 500m for the line coordinate and every 80m for the receiver coordinate. In the complex zones these spacings were smaller. The next step is to calculate the weights that have to be specified for the objective function used for reflection tomography. Our weights are a function of distance between the selected shot positions and receiver positions and the estimated residue associated to each set of traveltimes from one reflector. The latter one is obtained by taking the mean value for the difference between the set of picked traveltimes and the set of ideal CMP hyperbolas as calculated for the CMP difference plot.

The first set of traveltimes

The final set of traveltimes used for the first iteration of the SMART method was constituted as follows: For each reflector a list of selected traveltimes was established as function of the regularity of the CMP difference plot. For each set an estimation of the residue was made and the associated weights were calculated. At the end of this phase the traveltimes of 12 reflections were picked. The selected events are depicted in Figure 1. Those indicated with a dashed line were difficult to pick.

THE SMART ITERATIONS

After completion of the preparative phase, the phase in which the ultimate accuracy for both velocity and structural model is sought can be started. Our strategy is based on the SMART method. This interpretive prestack approach requires the following steps:

a Processing of the traveltimes by reflection tomography l Depth migration of the seismic data or application of the PICLI method (See Part 1 of this paper) with the model obtained in the previous step l a a l Interpretation of the prestack migrated data. Try to find important geological events which aren't flat in the coherency panels. Computation of the traveltimes associated with the picked events in the migrated data by means of raytracing. Adding these traveltimes to the already existing set of traveltimes Repeat this sequence until a satisfactory depth image is obtained; hence the name of the method Sequential Migration Aided Reflection Tomography. In complex structures like this one a number of iterations is needed in order to find the minimum set of traveltimes that allows for an accurate determination of the velocity model.

We executed a number of iterations. Each iteration improved the velocity model as could be concluded from the increasing flatness of the coherency panels and the more and more geological image that was obtained after stacking the prestack depth migrated common offset gathers.

First iteration

For first iteration only the most reliable traveltimes as picked in the preparative phase were used. For this reason the traveltimes of the deep 'base of the salt' reflector were not included. The picking of these traveltimes in the unmigrated data was quite difficult and it was expected that after the first iteration of the SMART method they could be better determined in the cube of migrated common offset gathers. As a consequence the velocity model determined by tomography was only correct for the regions above the salt. However the effect of the migration on the deeper interfaces allowed us to pick in the cube the events associated with the base of the salt. After obtaining the associated traveltimes, it was possible to add to the set of earlier picked traveltimes, a great number of traveltimes associated with the base of the salt.

Second iteration

With these additional traveltimes -departing from the previous modela new velocity model was calculated. The velocity model is depicted in Fig 4 . The lines are those parts of the reflectors that are determined by the tomography. The data were migrated. The post migration stack is shown in Figure 5. For the upper part, as could be expected the result of migration didn't change much. However the base of the salt became now more visible and at the left side where a maximum of traveltimes could be furnished, a very good definition of the subsurface was found. The mean residues for most of the inverted reflectors are very low (2 msec to 10 msec). In regions where no or little traveltimes could be prepared the model and thus the result of migration is rather disappointing. However when comparing this result with the conventional stack, it can be easily concluded that even in this case, picking of the base of the salt as a continuous event is easier than in any time section. The coherency panels are presented in Figure 6. The panel at 2000m indicates that the velocity model is quite correct. The event D is slightly dipping, probably because of the little number of traveltimes that are available at the left border of the dataset. The event G is flat. The strongly dipping events in between them are multiples. Deeper there are a number of more or less flat events, that correspond with the base of the salt. In the central part of the data, where no traveltimes were available, the model is highly undetermined and as a consequence the coherency panel at 8000m doesn't show any flat event. Between 3000 and 4000m depth some events are visible. They correspond with the base of the salt. Clearly there will be no image of this part of the base of the salt in the stack. At the right side of the dome at 10000 m the panel shows some continuous but dipping events. In that zone too little traveltimes were available for correctly determining the velocity model. The event marked P is one of the reflectors of the sedimentary sequence right of the dome. Somewhere below the base of the salt should be visible. It can be probably found around 3800m.

Interesting are the residues associated to the interfaces J and G. The residues on G are very low, corresponding with the flat events in the coherency panels. However at the interface J close to the dome the residues are relatively high. Therefore the traveltimes were corrected for this interface. The event has been repicked in the depth domain and better traveltimes were obtained after raytracing. Apart from a number of corrections, new traveltimes were added.

Third iteration

After the addition of these traveltimes it was expected to arrive at a very good definition of velocity model. It turned out that again the residues for almost all inverted times were very low. The model is shown in Figure 7. After migration of the data the individual common offset gathers were of a much better quality than those as obtained in the second iteration. However the data due to the presence of groundroll are rather noisy. This migrated noise has the tendency to mask in the central part the weak reflectors of the base of the salt. Therefore we decided to apply instead of a migration the PICLI method on these data. The stack together with a number of coherency panels is shown in Figures 8 and9. When comparing these results with those obtained in the second iteration, it is easy to observe the improvement of the imaging of the deeper part of the data. The salt base is correctly (flat and continuous) defined at the right side. This is especially due to the additional traveltimes associated with the sediments right of the dome. The PICLI method removed migration artefacts between 6000 and 9000m which made the continuation still more pronounced. At the left side where we added intermediate traveltimes associated with the unconformity the presence of a small fault in the basement can be observed.

In the coherency panel at 2000m there is now a clearly continuous event. Above the salt event there aren't so many changes, because we didn't change the traveltimes neither the a priori information. The most striking difference can be found in the coherency panel at 8000m. We improved the traveltime information of the upper sediments and we added traveltimes left and right of the dome. As a consequence great parts of the model between 6000 and 9000m are now better determined yielding a better migrated image. The panel at 10000m gives also a better image of the base of the salt. There are not so much improvements for the reflectors between 1200 and 2200m. This is probably related to extend of the corrections applied on the traveltimes in the second iteration.

The interpretability of the migrated data has become excellent, so that in a next iteration it will be possible to add still a few more traveltimes for those reflectors for which the coherency panels aren't completely flat.

CONCLUSIONS

Using a 2D dataset covering a really complex structure we have shown that standard processing cannot resolve the imaging problems related to complex wave propagation in such situations.

We have demonstrated the feasibility of the velocity determination method called SMART when applied on real data.

ACKNOWLEDGMENTS

The data for this article were kindly provided by Shell. This research was rried out as part of the Prestack Structural Interpretation Consortium Project (PSI). The authors hereby acknowledge the support provided by e sponsors of this project.

  Sinoquet, and Bertrand Duquet, Institut Francais du Pe'trole, France

Figure 1

 1 Figure 1 Unmigrated partial stack. The events for which traveltimes are picked are Indicated by solid or dashed lines according to the reliability of the picking.

Figure 2 AFigure 3

 23 Figure 2 A part of a CMP gather at 6000m Note the non-hyperbolic character of the reflector al 2300 msec.

4

The SMART method: Application


domain: math sdu

text:
Fabrice Jurado

Delphine Sinoquet

Andreas Ehinger

3D reflection tomography designed for complex structures

S u m m a r y

A 3D reflection tomography that can determine correct subsurface velocity structures is of strategic importance for an effective use of 3D prestack depth migration. We have developed a robust and fast 3D reflection tomography that is designed to handle complex models. We use a B-spline representation for interface geometries and for the lateral velocity distribution within a layer and we restrict the vertical velocity variation to have a constant gradient. We solve the ray tracing problem by use of a bending method with a circular ray approximation within layers. For the inversion we use a regularized formulation of reflection tomography which penalizes the roughness of the model. The optimization is based on a quadratic programming formulation and constraints on the model are treated by the augmented Lagrangian technique. We show results of ray tracing and inversion on a rather complex synthetic model.

Introduction

 Ehinger and Lailly, 1995, haveshown the interest of reflection tomography for computing velocity models adequate for the seismic imaging of complex geologic structures. In 2D, reflection tomography has proved its effectiveness in this context [START_REF] Jacobs | Velocity model determination by the SMART method[END_REF]. In 3D, Guiziou et al.,199 1, have developed a ray tracing based on a straight line ray approximation within a layer and an inversion of poststack data. But it suffers of a non derivability of its traveltime formula due to the Gocad interface representation. We describe a 3D tomography that handles models with the necessary derivability and allows inversion of more complex kinematics by the use of a more accurate traveltime calculation.

Model description

We choose a blocky model representation of the subsurface, each layer being associated with a geological macrosequence. A velocity law has to be associated with each layer (Figure 1). The form of the velocity law is = y) + where is the lateral velocity distribution (described by cubic B-spline functions) and k is the vertical velocity gradient. Using blocky models can lead to difficulties associated with the possible non-definition of the forward problem (situations where there is no ray joining a source to a receiver) and more generally to all kind of difficulties involved in discontinuous kinematics. The blocky model representation allows velocity discontinuities as they exist in the earth and thus to straightforwardly integrate a priori information on velocities (see [START_REF] Lailly | Smooth velocity models in reflection tomography for imaging complex structures[END_REF], for a general discussion on blocky versus smooth model s for seismic imaging of complex geologic structures). We use a cubic B-spline representation for interface geometries. This Fig. 1 Example of a model presenting a pinch-out. Continuous lines are the reflecting parts of interfaces, the dashed line is a virtual part that has no physical meaning. Specifying a ray signature allows to sort good and bad rays.

gives the (at least) regularity required for proper ray tracing calculations. Following [START_REF] Clarke | Ray shooting in 2D and 3D complex, blocky media[END_REF], the B-spline interfaces are defined in the whole domain and can cross each other. The part of an interface where the velocity is the same on both sides is called a virtual part (Figure 1); this part exists only for convenience, it has thus no physical meaning.

Ray tracing method

We compute normal (zero-offset) rays as well as multi-offset rays. We have to define the ray paths in which we are interested in. This consists in specifying the ray signature which is the ordered sequence of interfaces crossed by the ray before and after reflection. In other words, specifying the ray signature means specifying the nature of the seismic event we are interested in. In particular it can be used to remove artificially reflected rays i.e. rays that are reflected on the virtual part of an interface (Figure 1).

Two-point ray tracing

The two-point ray tracing method (zero or multi-offset) is inspired from the classical bending method. It consists of: -a ray initialization, i.e. we choose a trajectory that starts at the source, hits the reflector and ends at the receiver satisfying a given signature. This trajectory is defined by the successive impact points on the crossed interfaces; a formula for calculating the traveltime between two successive impact points; -the bending itself, i.e. the optimization technique that moves impact points until the total traveltime is stationary (Fermat's principle).

Ray initialization

For a zero-offset ray, the initial trajectory for the bending method is obtained by finding the successive intersection points of all interfaces with a circle (why a circle? see below...) starting from the source, with a given tangent vector (i.e. shooting angle), so as to reach the nearest interface. We apply Snell's law at the impact point to simulate the transmission at the interface, and iterate the circleinterface intersection procedure until the reflector is reached. The ray initialization with a given tangent vector can be very useful in the case of multiple arrivals in order to guide the bending method so as to find the desired arrival. If the tangent vector direction is not specified then by default the vertical direction is used. In the case of multi-offset CMP acquisition, an offset continuation strategy is applied. A new ray is initialized by copying, from the previously calculated ray, the sequence of impact point coordinates except for the source and receiver locations which are updated according to the considered offset.

Traveltime formula

For the calculation of the traveltime between two successive impact points, [START_REF] Guiziou | 3D Ray tracing through complex triangulated surfaces[END_REF], have proposed a formula based on a straight line ray approximation. To improve the accuracy of the traveltime calculation for velocity distribution with constant vertical gradient (see ray tracing results, Figure 5), we use a formula suited for such a velocity. In the case of constant lateral velocity, it is well known that the ray trajectory is a circle that lies in the vertical plane defined by the two successive impact points and P n . The traveltime t n is calculated by (see [START_REF] Virieux | Propagation in inhomogeneous media, Ray theory[END_REF]:

where angles and are explained on Figure 2. Note that we have translated the coordinate system origin to P n-1 and rotated the (x,y) plane to be in the vertical plane in which the ray lies (i.e. we deal with a 2D problem). Then the circle center depth is = and the coordinate is obtained by finding the intersection between this plane and the median of segment In the case of lateral velocity variations, is approximated by the mean of the velocity values at and M (Figure 2). In the case of strong lateral velocity variations, this approximation can be poor. In order to improve the accuracy of the computed traveltime, we can then introduce ghost interfaces within the layer. The ghost interface has no physical meaning: its role is just to allow an extra impact point, i.e. one more degree of freedom, so as to make the approximation more accurate.

Bending method

We define the total traveltime between the source and the receiver b y = which is a differentiable function of impact point coordinates (thanks to

We look for impact point locations that make the total traveltime minimum. This is done using a quasi-Newton method which requires only calculation of the first derivative of the traveltime with respect to impact point coordinates. This method has asymptotically the same convergence rate as a Newton method but is much cheaper and also much more stable.

Since interfaces can cross each other, the ray signature may change from one iteration of the bending to another (whenever an impact point moves across another interface) and the bending algorithm can find rays that violate the specified signature. We allow the bending to find such rays but we discard them by checking the signature a posteriori. Thus the bending algorithm is not slowed down by extensive signature management and at the same time we are sure to discard erratic rays. Note that finding no ray for a given sourcereceiver pair is quite acceptable since we cannot be sure that there exists such ray when dealing with a blocky model. Also, as in all bending methods, there is no guarantee, in the case of multiple arrivals, to find the arrival of interest: the solution depends on the ray initialization, An a posteriori check of the shooting angle of the ray allows to detect a change of the traveltime branch. Despite this minor limitation our ray tracing has been designed to be fast (thanks to the quite simple traveltime formula) and reliable (thanks to the B-spline representation of surfaces leading to the differentiable definition of the total traveltime). The bending algorithm converges towards a trajectory with zero traveltime derivatives with respect to impact point positions. This is an important condition to be fulfilled for correct calculation of the Jacobian.

Tomographic inversion method

The tomographic inverse problem consists in finding the modelnamely, the B-spline parameters of the lateral velocity and the interface depths for each layer (the constant velocity gradient has a given and fixed value) -that yields traveltimes that match the data. The data of the inversion are the prestack traveltimes. These data could be obtained for instance from zero-offset (poststack) traveltimes together with offset move-out information. We are able to invert multivalued traveltimes if different branches of a traveltime curve are defined by different shooting angles of the zero-offset ray. This a priori information can be read from the zero-offset data (first derivatives of the traveltime surface). Inversion of prestack traveltimes gives the possibility to invert more complex kinematics as compared to the direct inversion of stacking velocities. We also expect better stability since prestack traveltimes yield a bigger amount of information than the associated stacking velocity. As shown in Delprat-Jannaud and Lailly, 1993, the inverse problem has to be regularized. Indeed, some parts of the model may not be illuminated by rays and thus are not determined by the data. We represents the inverse of the uncertainty associated with traveltime Choosing the weights correctly comes down to finding the appropriate compromise between possibly contradictory pieces of information (see [START_REF] Lailly | Smooth velocity models in reflection tomography for imaging complex structures[END_REF]. The minimization of the non quadratic objective function is carried out with the Gauss-Newton method. In fact, for better stability, it is important to put constraints on the model. Each iteration thus consists in solving a constrained least-squares problem. We make use of a quadratic programming formulation which allows the inversion of a large number of data: with this formulation we just need to store the Hessian matrix instead of the Jacobian matrix. To handle equality and inequality constraints we make use of a specific implementation of the augmented Lagrangian method. We can thus, for example, fix the velocity at the surface or at well locations, or constrain the depth or the slope of an interface, or constrain the relative position of two interfaces.

Applications of ray tracing and inversion Synthetic model

We have created a synthetic model to test our ray tracing and tomographic inversion, The model is inspired by a real salt dome structure Fig. 3 Ray tracing on the synthetic model The depth of the ternary base vanes from 600m to 1800m and the depth of the chalk base varies from 600m to 2800m (see Figures 6 and7).

in the North Sea and consists in total of 7 interfaces. The upper two interfaces (Figure 3) are the tertiary base and the chalk base. Each interface is explicit in the depth direction (Z(x, y)) and represented by 30 x 15 B-spline coefficients on the independent axes. The model extends from -5000m to 13000m in x direction and from 0m to 9400m in y direction. The lateral variations of the velocity in the first layer are shown in Figure 6 (bottom right). The constant vertical velocity gradient is For the second layer, the lateral velocity is represented in Figure 7 (bottom right). The vertical gradient is Each lateral velocity distribution is represented by 14 x 8 B-spline coefficients. The acquisition survey consists in zero-offset shots and CMP gathers. 360 x 188 zero-offset shots are positioned on a regular grid of 50m width. We place 37 x 24 CMPs on a regular grid of 400m width with 60 offsets varying from 50m to 3000m The zero-offset shots were first calculated by an accurate ray tracing developed by Clarke et al., 1994, in order to access the shooting angles in a simple way rather than by an interpretation. Using these shooting angles we then recomputed the zero-offset data and the multi-offset data with our ray tracing.

Ray tracing results

In Figure 5, we present the improvement of the circular ray formula with respect to a straight line ray approximation. Figure 5 (left) shows results of ray tracing using the circular ray approximation compared to a straight line ray approximation, i.e. in each layer the ray is a line and the traveltime is computed using the mean slowness value (see [START_REF] Guiziou | 3D Ray tracing through complex triangulated surfaces[END_REF]. We see that the difference can be large especially for large offsets. In fact, the difference would be larger when dealing with more complex interfaces. The ray tracing with the circular ray approximation is faster than that with the straight line ray approximation when no ghost interfaces are used; moreover in practice one ghost interface is necessary for the straight line ray approximation to deal with vertical velocity gradient and increase the accuracy, its CPU time is then at least multiplied by 3. To certify our ray calculation, Figure 5 (right) compares results between the ray tracing with circular ray approximation and the ray tracing of [START_REF] Clarke | Ray shooting in 2D and 3D complex, blocky media[END_REF]. We see the very good agreement even for large offsets. For this common mid-point, the exact ray tracing found three different branches but two of them are very close, so we have fed the bending method with only two really different shooting directions and have found only two branches. On Figure 4, the results of zero-offset raytracing with the shooting method and the bending method (circular ray approximation) show also that the bending method can have difficulties to find the desired arrivals when the initialization shooting directions are too close for two branches. The advantage of our ray tracing is to be 20 times faster (zero-offset shots) than the exact ray tracing.

3D reflection tomography

Fig. 4 Zero-offset times associated with the chalk base of the synthetic model along a line y = 3025m They were obtained with a shooting method [START_REF] Clarke | Ray shooting in 2D and 3D complex, blocky media[END_REF] and with the bending method (with the circular ray approximation). The bending method retrieved two branches of the traveltime curve but was not able to retrieve the third branch due to an initialization direction too close to that of the second branch which then was calculated twice.

Inversion results

We performed inversion of the zero-offset and multi-offset traveltime data associated with the two first layers of the synthetic model by a layer-stripping approach. First, we invert for the tertiary base starting from an initial constant horizontal interface and an initial constant velocity; it leads to an inversion for 450 B-spline coefficients for the interface and 112 coefficients for the velocity. The model was regularized by means of second derivatives without any a priori model. Also weights on the traveltimes data were tuned to give sufficient importance of the large offset data with respect to the (huge) number of zero (or small) offset data. We invert almost 130000 traveltime data.

Figure 6 shows the exact tertiary base compared with the interface obtained after 17 iterations of the inversion and the same comparison between velocities. The RMS traveltime misfit is 0.2ms (only 5 iterations were necessary to obtain a RMS misfit smaller than lms). We now fix the tertiary layer to the obtained result and invert the chalk base and the lateral velocity, again starting from constant values. Figure 7 shows the exact chalk base compared with the interface obtained after 14 iterations and the same comparison between velocities. The RMS traveltime misfit is 0.7ms. Only a few traveltime misfits remain large: they correspond to a change of branch of the traveltime curve, the initialization shooting angles being not sufficient to distinguish the different arrivals (see Fig, 5 Triplicated multi-offset times associated with the chalk base of the synthetic model for a CMP ( x = 6775m) y = 4125m). Left side: comparison between our bending method with the circular ray approximation and a bending method with a straight line ray approximation. Right side: comparison between our bending method with the circular ray approximation and a shooting method [START_REF] Clarke | Ray shooting in 2D and 3D complex, blocky media[END_REF]. Two initialization angles for the bending provide the two branches.

section on ray tracing results).

The CPU time amounts to 20 and 45 minutes for one iteration of the inversion for layer 1 and 2 (respectively) on a Silicon Graphics Indigo² workstation (including ray tracing).

Conclusions

The 3D reflection tomography described in this paper has proven to be efficient for a geologic synthetic structure. Using B-spline representation of the interfaces and of the velocity gives the necessary regularity for the ray tracing and the inversion. The ray tracing relies on a bending method which is able to compute multiple arrivals providing that we initialize the process with reliable shooting angles (first derivatives of zero-offset data). The bending method combined with the circular ray approximation allows fast calculation and gives a sufficient accuracy even with large velocity variations. The regularized inversion formulation with constraints is robust and can overcome the difficulties associated with irregular model illumination by rays which is particularly true for complex structures.

Fig. 2

 2 Fig. 2 Description of the traveltime calculation for circular ray path. C is the center of the circle, and are successive impact points. is the mean of the velocity values at , k is the vertical velocity gradient.

  have chosen to regularize by penalizing the second derivatives of the difference between the actual model and an a priori model. We thus solve the following minimization problem. where NT is the number of traveltime data and the model m is composed of NV velocity laws and of NZ interface depths The first term in the objective function measures the misfits between the traveltimes computed by the ray tracing in the current model m and the traveltime data. The regularization term means that we are looking for a model with a roughness close to the a priori model or, if there is no a priori model, a small roughness. represent weights which calibrate in the objective function the different terms according to their physical meaning.

Acknowledgements

This research was carry out as part of the Kinematic Inversion Methods consortium project (KIM). The authors hereby acknowledge the support provided by the sponsors of this project. We thank Elf

Interface result

Aquitaine for providing us with models used to create the synthetic model. We thank Richard Clarke for the interesting comparisons done with his ray tracing, Fatmir Hoxha for interesting discussions on circular ray calculation and Patrick Lailly for his fruitful comments.

Synthetic interface


Reflection tomography is an efficient method to determine a subsurface velocity model that best fits the traveltime data associated to the main events picked on the seismic sections. A careful choice of the model representation has to be done: a blocky model representation based on regular gridded bspline functions has been proposed. This flexible parameterization allows accurate and robust inversion but can lead to a huge number of parameters. An adaptive parameterization that enables to account for local complexities and inhomogeneous ray coverage is considered.

Introduction

Reflection tomography allows to determine a subsurface velocity model from the traveltimes associated with the main reflections of the seismic waves. Solving the forward and inverse problems requires a careful choice of the model parameterization: it affects the computation cost for the ray tracing and has strong consequences on the robustness of the inversion. The KIM consortium proposed a blocky model representation [START_REF] Clarke | Raytracing and tomography friendly geological models : KIM 1996 Annual Report[END_REF] based on b-spline functions: its effectiveness has been proven on several applications on complex data [START_REF] Ehinger | 3D tomographic velocity model determination for two North Sea case studies : 63 rd Conference and Technical Exhibition[END_REF]. Nevertheless, the use of regular grids to define b-spline functions leads often to the inversion of useless parameters, especially in case of models with local complexities or with inhomogeneous ray coverage. Several authors have developed adaptive model parameterizations for tomography: [START_REF] Böhm | 3D adaptive tomography using Delaunay triangles and Voronoi polygons[END_REF] and [START_REF] Cox | Data-driven tomographic inversion of focusing operators : 71 st[END_REF] have proposed automatic data-driven adjustments of Delaunay triangle parameterization. The purpose of our work is to provide an adaptive model parameterization, based on b-spline functions, in order to limit the number of inverted parameters and consequently to make the inversion more efficient.

Formulation of reflection tomography

The forward problem of our reflection tomography [START_REF] Jurado | 3D reflection tomography designed for complex structures : 66 th Ann[END_REF] is a two-point ray tracing based on a bending method: complex kinematics (converted waves, multi-valued arrivals, multiple reflections but also direct arrivals) associated with arbitrary acquisition surveys (classic surface data, OBC, VSP, walk away) can be modeled. The inverse problem consists in the minimization of a least-square objective function is a regularization term, composed of the sum of second derivatives of the interface depths and velocity variations. This penalization term is necessary to stabilize the inversion and allows the introduction of a priori information.

In our approach, the subsurface velocity is represented by a blocky model with interfaces described by 2D explicit b-spline functions and layer velocities described by 2D or 3D explicit b-spline functions. Constant anisotropy parameters (assuming transverse isotropy) in each layer can also be inverted (see [START_REF] Broto | Anisotropic traveltime tomography for depth consistent imaging of PP and PS data[END_REF]. The choice of explicit functions and regular discretization grids is crucial for fast ray-tracing calculations. In order to easily handle the model changes during the inversion process, we define the interfaces and the velocities on the whole definition domain of the model, creating then virtual parts of interfaces or velocities (virtual meaning non-physical). Note that this parameterization can be used to represent hybrid models composed of velocity blocks delimited by reflectors, and of smoothly varying velocities with embedded reflectors. Examples of complex solution models found by tomography applied on real data are displayed in Figure 1.

A realistic 3D model may require up to 50000 b-spline parameters. [START_REF] Delbos | Trust-region Gauss-Newton method for reflection tomography[END_REF] have developed sophisticated optimization techniques for the inversion of big models which allow the introduction of linear constraints on the model. An other way to improve the efficiency of the inversion of such models is to work on the parameterization itself to limit the number of parameters.

The parameterization described above requires the a priori choice of a discretization step of the bspline functions: for instance, the ideal discretization step of the interface represented on figure 2 would be small in the complex zones and large in the simple areas (figure 2 right). Such an adaptive parameterization allows to drastically reduce the number of parameters used to describe the interface (compared to a regular grid, figure 2 left). More generally, the approach we propose here aims at reducing the inverted model size by removing from the inversion the useless, that is to say the nonphysical (e.g. virtual parts of reflectors), non-illuminated, or redundant b-spline parameters to improve the stability of the inversion, and to reduce both memory use and computation costs.

Macro-parameter methodology

To reach these objectives, we introduce macro-parameters, which are linear combinations of model parameters. Defining macro-parameters results in replacing many useless parameters (non-physical, non-illuminated, or redundant information in simple areas) by few macro-parameters to obtain a reduced sized model. In this way, an irregular parameterization can be obtained from the dense discretization (figure 2, left) by replacing one or several groups of parameters in basic macroparameters as constant values or planes (i.e. constant or linear extrapolation to deal with virtual parts of the model or unilluminated areas), or in more sophisticated ones as the merging of these parameters in one of their embedded b-spline parameterization (to deal with local complexities).

The choice whether to replace a group of parameters by a macro-parameter or not, must be carried out periodically in the inversion procedure to take into account model variations. It can be done automatically through the computation of update indicators. These indicators estimate the coarser bspline discretization that explains the data complexity, and are classically based on the resolution matrix analysis, or some related approaches, which are often expensive and a posteriori methods [START_REF] Böhm | 3D adaptive tomography using Delaunay triangles and Voronoi polygons[END_REF][START_REF] Cox | Data-driven tomographic inversion of focusing operators : 71 st[END_REF][START_REF] Cox | Data-driven tomographic inversion of focusing operators : 71 st[END_REF]. We are currently developing a cheaper indicator (based on the gradients), that gives the optimal a priori b-spline discretization to fit the data (as illustrated in the next sections).

Application to realistic synthetic data

The method described in the previous section is here applied on a realistic synthetic model composed of 8 interfaces and 3 velocities, including a smooth (3D) velocity located under the second interface (figure 4). The six deepest interfaces (3024 parameters) and the 3D velocity (6000 parameters) are inverted simultaneously with the classic and macro-parameter version of the inversion software. The deepest interfaces are only partly illuminated. We replace the many b-spline parameters in the nonilluminated zones (hatched zones on the figure 3) by macro-parameters describing simple planes. Besides, the 3D velocity is defined all over the model, even if it only applies to the volume located under the second interface, and is not illuminated beneath the deepest interface (4 km). We thus replace the b-spline parameters located above 1 km and below 4 km by constant macro-parameters.

The inversion with macro-parameters furnishes (figure 4) very satisfactory results (RMS = 0.8 ms), with a reduced number of model parameters (figure 5). These results are equivalent to the one obtained with the classic version, but with significant benefits: the macro-parameter version requires less memory, and is 25% faster. Moreover, the tuning of the regularization weights is considerably simplified, because of the introduction, via macro-parameters, of hard constraints in the areas of poor ray coverage, which may be seen as an induced (model intrinsic) regularization.

Application to a 2D model with local complexities

To illustrate the effectiveness of our update indicator we now apply our method to synthetic data calculated on a model composed of a velocity and a locally complex interface. We compute the gradient-based indicator for the interface of the initial model (plane interface at z = 800 m and constant velocity) in order to determine the b-spline discretization step required to fit the data (figure 6). A discretization step is eligible if its indicator is near zero. We notice that the interface is locally complex especially from 6 to 10 kilometers where a discretization step of 78 meters is required. The classic version is consequently performed with this step, leading to the inversion of 131 b-spline parameters. However, between 0 and 6 kilometers, which is relatively simple region, it is possible to choose a discretization step of 156 meters (312 meters is not good enough). We have run the macroparameter version to invert the data with the chosen discretization step: 156 meters between 0 and 6 kilometers and 78 meters between 6 and 10 kilometers. The inversion is performed with 95 parameters instead of 131. The solution model (figure 7) obtained with this parameterization is equivalent to the one obtained with the classic version but with a significant gain (27%) on the number of inverted parameters.

Conclusions

We have developed a new adaptive parameterization approach, based on an original formalism, that makes our reflection tomography even more efficient. Our model description now combines the efficiency of regular parameterization for ray-tracing and the flexibility of irregular grids for inversion.

The results obtained with the macro-parameter approach compared to the classic version of the software show multiple benefits: the inversion is more stable and less sensitive to regularization weights settings (this is mainly due to the model intrinsic regularization), the computing costs and memory use are significantly reduced. This method is mainly automatic thanks to the computation of appropriate update indicators. Then, macro-parameters coupled with adequate update indicators allows an adaptive parameterization for reflection tomography. This formalism may also be used for the introduction of hard a priori constraints on the model for underdetermined problems. Iteration number

Figure 2 :Figure 3 :

 23 Figure 2: Regular versus irregular parameterization of an interface with local complexities. Many parameters on the left are useless and are not worth to be inverted.

Figure 4 :

 4 Figure 4: Solution model obtained by the macroparameter version. The bold lines correspond to the illuminated parts of the interfaces.

Figure 5 :

 5 Figure 5: Inverted model size for the macroparameter (circles) and classic (squares) versions as a function of tomography iterations.

Figure 6 :

 6 Figure 6: Update indicator. The discretization step that best fits the data gives a near-zero indicator.

Figure 7 :

 7 Figure 7: Solution model obtained by the macroparameter version with two different b-spline discretization steps.

Z-99 Flexible b-spline model parameterization designed for reflection tomography

T. PERDRIZET , D. SINOQUET

Institut Français du Pétrole, 1&4 avenue de Bois-Préau 92852 Rueil-Malmaison Cedex -France