Method of estimating elastic and compositional parameters from seismic and echo-acoustic data

• US 6,876,928 B2
• Filed: 09/26/2003
• Issued: 04/05/2005
• Est. Priority Date: 10/22/1999
• Status: Expired due to Fees

• Current Assignee: CGG Services (NL) BV (CGG)
• Original Assignee: Jason Geosystems B.V. (Fugro N.V.)
• Inventors: Debeye, Hendrik Willem Johan, Van Riel, Paul
• Primary Examiner(s): Hoff, Marc S.
• Assistant Examiner(s): Gutierrez, Anthony
• Application Number:
  US10/672,751
• Publication Number:
  US 20040064294A1
• Time in Patent Office: 557 Days
• Field of Search: 702/2, 702/16, 702/17, 702/189
• US Class Current: 702/2
• CPC Class Codes: G01V 1/306 for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles

Abstract

A method for determining from measured reflection data on a plurality of trace positions, a plurality of subsurface parameters. The method includes the steps of: preprocessing the measured reflection data into a plurality of partial or full stacks; specifying one or more initial subsurface parameters defining an initial subsurface model; specifying a wavelet or wavelet field for each of the partial or full stacks of the measured reflection data; calculating synthetic reflection data based on the specified wavelets and the initial subsurface parameters; optimizing an objective function, including the weighted difference between measured reflection data and synthetic reflection data for a plurality of trace positions simultaneously; and outputting the optimized subsurface parameters. A device for implementing this method is also included.

Claims

1. A method for determining from measured reflection data on a plurality of trace positions one or more subsurface parameters, said method comprising the steps of:

(a) preprocessing the measured reflection data into a plurality of partial or full stacks;

(b) specifying one or more initial subsurface parameters defining an initial subsurface model;

(c) specifying a wavelet or wavelet field for each of the partial or full stacks of the measured reflection data;

(d) calculating synthetic reflection data based on the specified wavelets and the initial subsurface parameters;

(e) optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data for a plurality of trace positions simultaneously; and

(f) outputting the optimized one or more subsurface parameters.

2. The method according to claim 1, wherein the step of optimizing the objective function comprises minimizing the objective function:

$F_{\mathrm{residual}}=\underset{i=1}{\overset{\#<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{stacks}}{\sum <br><br>}}<br><br>w_{\mathrm{residual},i}<br><br>\underset{j=1}{\overset{\#<br><br>\mathrm{traces}}{\sum <br><br>}}<br><br>L_{P,\mathrm{residual}}<br><br>\left(S_{i,j}-m_{i,j}\right)$

wherein s_i,j represents trace j of measured reflection data for stack i, m_i,j represents trace j of the modeled synthetic reflection data for stack i, w_residual,i represents a weighting factor for stack i, #traces represents the total number of traces, #stacks represents the total number of stacks and L_P,residual is an adjustable norm of the residuals (s_i,j−

m_i,j), which provide the deviation of the synthetic reflection data from the measured reflection data.

3. The method according to claim 2, wherein the objective function comprises one or more stabilization terms and/or one or more correction terms.

4. The method according to claim 3, wherein a stabilization term is a measure for the deviation of the reflectivity away from 0.

5. The method according to claim 4, wherein said measure comprises $F_{\mathrm{reflectivity}}=\sum _{}^{}$

i=1#⁢

 

⁢

stacks⁢

wreflectivity,i⁢

∑

j=1#⁢

traces⁢

LP,reflectivity⁡

(ri,j)

wherein r_i,j is the reflectivity trace for stack i and trace j, w_{reflectivity,i} is a weighting factor for stack i, #traces is the total number of traces, #stacks is the total number of stacks and L_{P,reflectivity} is an adjustable norm of the reflectivities.

6. The method according to claim 3, wherein a stabilization term is a measure for the parameter contrast.

7. The method according to claim 6, wherein said measure comprises $F_{\mathrm{contrast}}=\sum _{}^{}$

n=1#⁢

 

⁢

parameters⁢

wcontrast,n⁢

∑

j=1#⁢

traces⁢

LP,contrast⁡

(cj,n)

wherein c_j,n is the contrast trace for the n^th subsurface parameter, w_contrast,n is a weighting factor for the n^th parameter, #traces is the total number of traces, #parameters is the number of parameters, and L_P,contrast is an adjustable norm of the contrasts.

8. The method according to claim 7, wherein a stabilization term is a measure for the deviation of the calculated subsurface parameters from a priori specified functional relationships between subsurface parameters.

9. The method according to claim 8, wherein said measure comprises $F_{\mathrm{functions}}=\sum _{}^{}$

v=1#⁢

 

⁢

functions⁢

wfunctions,v⁢

∑

j=1#⁢

 

⁢

traces⁢

LP,functions⁢

fv(pj,…

⁢

 

,⁢

pj,#⁢

 

⁢

parameters)

wherein f_v represents the deviations of the subsurface parameters at trace j away from the v^th functional relationship between different subsurface parameters, w_functions,v is a weighting function for the v^th functional relationship, #traces is the number of traces, #functions is the number of functional relations and L_P,functions is an adjustable norm of said deviations.

10. The method according to claim 3, wherein a stabilization term is a measure for the deviation of the subsurface parameters from the initial subsurface parameters.

11. The method according to claim 10, wherein said measure comprises $F_{\mathrm{initial}}=\sum _{}^{}$

n=1#⁢

 

⁢

parameters⁢

winitial,n⁢

∑

j=1#⁢

traces⁢

LP,initial⁡

(pj,n-pinitial,j,n)

wherein (p_j,n−

p_initial,j,n) represents the trace with the difference between the calculated subsurface parameter n at trace position j and the corresponding initial subsurface parameter, w_initial,n is a weighting factor for the n^th parameter, #traces is the total number of traces, #parameters is the number of parameters and L_P,initial is an adjustable norm of said difference.

12. The method according to claim 3, wherein a stabilization term is a measure for the lateral variability of the parameters.

13. The method according to claim 12, wherein said measure comprises $F_{\mathrm{lateral}}=\sum _{}^{}$

n=1#⁢

 

⁢

parameters⁢

∑

l=1#⁢

 

⁢

neighbors⁢

∑

j=1#⁢

traces⁢

wlateral,n⁡

(rj,l)⁢

LP,lateral⁡

(dj,l,n)

wherein d_j,n,l is the difference of the samples of parameter p_n at traces j and l, corrected with any difference in the initial model, w_lateral,n (τ

_j,l) is a trace for parameter n describing the weighting for each parameter sample where this weighting is a function of τ

_j,l which is a trace which at each parameter sample provides a measure of the local correlation between the traces j and l, #traces is the number of traces, #neighbors is the number of neighboring traces used in the calculation, #parameters is the number of parameters and L_p,lateral is the adjustable norm of said differences d_j,n,l.

14. The method according to claim 13, wherein the parameter difference d_j,l,n is defined as

(d_j,l,n)(t_k)=p_l,n(t_k+Δ

t_j,l,k)−

p_j,n(t_k)−

(p_initial,l,n(t_k+Δ

t_j,l,k)−

p_initial,j,n(t_k))

wherein Δ

t_j,l,k is the time shift at parameter sample k which time aligns the parameters of trace l to trace j at sample k, where surrounding trace samples are interpolated if at time t_k+Δ

t_j,k,l a sample is not defined, where r_j,l is now the local correlation incorporating the time shift and L_P,lateral is an adjustable norm on the parameter differences.

15. The method according to claim 2, wherein the adjustable norm L_P of a variable x for variable x comprises $L_p$

(x)=[∑

k=1#⁢

 

⁢

samples⁢

(

xk

)P]

wherein P is a user adjustable factor.

16. The method according to claim 15, wherein the norm L_P is normalized by or with at least one of:

(i) exponentiation with 1/P;

(ii) the number of samples; and

(iii) the square root of the variance of x.

17. The method according to claim 1, wherein a correction term is a measure for the differential time shifts between traces of measured reflection data stacks.

18. The method according to claim 17, wherein said measure comprises $F_{\mathrm{time}}=\sum _{}^{}$

i=2#⁢

 

⁢

stacks⁢

wtime,i⁢

∑

j=1#⁢

traces⁢

LP,time(τ

i,j-τ

 

0⁢

i,j)

wherein τ

_0,i,j is the trace with the initial time values of the time stretch and squeeze control points for stack i and trace j and τ

_i,j is the time of shifted control points, w_time,i is a weighting factor for stack i, #stacks is the number of stacks, #traces is the number of traces and L_P,time is an adjustable normalization factor of the difference between τ

_0,i,j and τ

_i,j.

19. The method according to claim 1, wherein a stabilization term is a measure for the parameter differences between reflection data acquisition surveys taken at different points in time.

20. The method according to claim 19, wherein said measure comprises $F_{\mathrm{timelapse}}=\sum _{}^{}$

k=2#⁢

 

⁢

surveys⁢

wsurvey,k⁢

∑

n=1#⁢

 

⁢

parameters⁢

wparameters,n⁢

∑

j=1#⁢

 

⁢

traces⁢

LP,timelapse(pj,n,k-pj,n,k-1)

wherein (p_j,n,k−

p<sub>j,n,k−

1</sub>) is the difference in parameters of trace j and parameter n between survey k and its time preceding survey k−

1, w_parameters,n is a weighting factor for each parameter n, w_survey,k is a weighting factor for each survey k, L_P,timelapse is an adjustable norm of said difference, #surveys represents the number of surveys and #parameters represents the number of parameters.

21. The method according to claim 1, wherein constraints and/or constraint functions are applied to one or more of the subsurface parameters;

and/or are applied to control changes of subsurface parameters between surveys taken on different times; and

/or constrain outside a specified subsurface zone the parameter changes between surveys to small values relative to changes expected within the specified subsurface zone; and

/or constrain the minimization by setting the subsurface parameters outside a specified subsurface zone from survey to survey at the same value.

22. The method according to claim 1, wherein for optimizing the objective function outside a specified subsurface zone only one or one set of different subsurface parameters are specified.

23. The method according to claim 1, comprising the generation of quality control information.

24. The method according to claim 23, wherein the quality information includes at least one of:

(i) synthetic data based on the optimized subsurface parameters;

(ii) residual data obtained by subtracting the synthetic data from the measured reflection data;

(iii) deviation data obtained by determining the deviations away from the initial subsurface parameters;

(iv) deviation data obtained by determining the deviations away from the corresponding functional relations; and

(v) deviation data obtained by determining the deviations away from well log data.

25. The method according to claim 1, wherein the reflection data is at least one of seismic data and time lapse data.

26. The method according to claim 25, wherein the time lapse data at each survey time comprises at least one seismic partial or full stack.

27. The method according to claim 25, wherein the seismic reflection data is determined from at least one of the following source-receiver combinations:

P-wave source and P-wave receiver, P-wave source and S-wave receiver, S-wave source and P-wave receiver, S-wave source and S-wave receiver.

28. The method according to claim 1, wherein the subsurface parameters comprise elastic parameters.

29. The method according to claim 28, wherein the elastic parameters comprise pressure wave velocities and/or shear wave velocities and/or densities in the subsurface and/or the subsurface parameters comprise any mathematical relation between pressure wave velocities and/or shear wave velocities, and/or densities.

30. The method according to claim 29, wherein the seismic data comprises at least two seismic partial or full stacks containing different angle dependant information on seismic reflections in the subsurface.

31. The method according to claim 1, wherein the subsurface parameters comprise compositional parameters representing the rock and fluid composition of the subsurface.

32. The method according to claim 1, wherein the reflection data is echo-acoustic data and the subsurface is human or mammal tissue or any other material.

33. A method for determining from measured reflection data on one or more trace positions a plurality of subsurface parameters, said method comprising the steps of:

(a) preprocessing the measured reflection data into a plurality of partial or full stacks;

(b) specifying a plurality of initial subsurface parameters defining an initial subsurface model;

(c) specifying a wavelet or wavelet field for each of the partial or full stacks of the measured reflection data;

(d) calculating synthetic reflection data based on the specified wavelets and the initial subsurface parameters;

(e) simultaneously optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data, for said plurality of subsurface parameters and for each trace position separately; and

(g) outputting said plurality of optimized subsurface parameters.

34. The method according to claim 33, wherein the step of optimizing the objective function comprises minimizing the objective function:

$F_{\mathrm{residual}}=\underset{i=1}{\overset{\#<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{stacks}}{\sum <br><br>}}<br><br>w_{\mathrm{residual},i}<br><br>\underset{j=1}{\overset{\#<br><br>\mathrm{traces}}{\sum <br><br>}}<br><br>L_{p,\mathrm{residual}}<br><br>\left(S_{i,j}-m_{i,j}\right)$

wherein s_i,j represents trace j of measured reflection data for stack i, m_i,j represents trace j of the modeled synthetic reflection data for stack i, w_residual,i represents a weighting factor for stack i, #traces represents the total number of traces, #stacks represents the total number of stacks and L_P,residual is an adjustable norm of the residuals (s_i,j−

m_i,j), which provide the deviation of the synthetic reflection data from the measured reflection data.

35. The method according to claim 34, wherein the adjustable norm L_P of a variable x for variable x comprises $L_p$

(x)=[∑

k=1#⁢

 

⁢

samples⁢

 

⁢

(

xk

)P]

wherein P is a user adjustable factor.

36. The method according to claim 35, wherein the norm L_P is normalized by or with at least one of:

(i) exponentiation with 1/P;

(ii) the number of samples; and

(iii) the square root of the variance of x.

37. The method according to claim 33, wherein the objective function comprises one or more stabilization terms and/or one or more correction terms.

38. The method according to claim 37, wherein a stabilization term is a measure for the deviation of the reflectivity away from 0.

39. The method according to claim 38, wherein said measure comprises $F_{\mathrm{reflectivity}}=\sum _{}^{}$

i=1#⁢

 

⁢

stacks⁢

wreflectivity,i⁢

∑

j=1#⁢

 

⁢

traces⁢

 

⁢

LP,reflectivity⁢

 

⁢

(ri,j)

wherein r_i,j is the reflectivity trace for stack i and trace j, w_{reflectivity,i} is a weighting factor for stack i, #traces is the total number of traces, #stacks is the total number of stacks and L_{P,reflectivity} is an adjustable norm of the reflectivities.

40. The method according to claim 37, wherein a stabilization term is a measure for the parameter contrast.

41. The method according to claim 40, wherein said measure comprises $F_{\mathrm{contrast}}=\sum _{}^{}$

n=1#⁢

 

⁢

parameters⁢

wcontrast,n⁢

∑

j=1#⁢

 

⁢

traces⁢

 

⁢

LP,contrastl⁢

 

⁢

(cj,n)

wherein c_j,n is the contrast trace for the n^th subsurface parameter, w_contrast,n is a weighting factor for the n^th parameter, #traces is the total number of traces, #parameters is the number of parameters, and L_P,contrast is an adjustable norm of the contrasts.

42. The method according to claim 41, wherein a stabilization term is a measure for the deviation of the calculated subsurface parameters from a priori specified functional relationships between subsurface parameters.

43. The method according to claim 42, wherein said measure comprises $F_{\mathrm{functions}}=\sum _{}^{}$

v=1#⁢

 

⁢

functions⁢

wfunctions,v⁢

∑

j=1#⁢

 

⁢

traces⁢

 

⁢

LP,functions⁢

 

⁢

fv⁢

 

⁢

(pj,1,⁢

 

⁢

…

⁢

 

⁢

pj,#⁢

 

⁢

parameters,)

wherein f_v represents the deviations of the subsurface parameters at trace j away from the v^th functional relationship between different subsurface parameters, w_functions,v is a weighting function for the v^th functional relationship, #traces is the number of traces, #functions is the number of functional relations and L_P,functions is an adjustable norm of said deviations.

44. The method according to claim 37, wherein a stabilization term is a measure for the deviation of the subsurface parameters from the initial subsurface parameters.

45. The method according to claim 44, wherein said measure comprises $F_{\mathrm{initial}}=\sum _{}^{}$

n=1#⁢

 

⁢

parameters⁢

winitial,n⁢

∑

j=1#⁢

 

⁢

traces⁢

 

⁢

LP,initial⁢

 

⁢

(pj,n-pinitial,j,n)⁢

 

wherein (p_j,n−

p_initial,j,n) represents the trace with the difference between the calculated subsurface parameter n at trace position j and the corresponding initial subsurface parameter, w_initial,n is a weighting factor for the n^th parameter, #traces is the total number of traces, #parameters is the number of parameters and L_P,initial is an adjustable norm of said difference.

46. The method according to claim 37, wherein a stabilization term is a measure for the lateral variability of the parameters.

47. The method according to claim 46, wherein said measure comprises $F_{\mathrm{lateral}}=\sum _{}^{}$

n=1#⁢

 

⁢

parameters⁢

∑

l=1#⁢

 

⁢

neighbors⁢

∑

j=1#⁢

 

⁢

traces⁢

wlateral,n⁡

(rj,l)⁢

LP,lateral⁡

(dj,l,n)

wherein d_j,n,l is the difference of the samples of parameter p_n at traces j and l, corrected with any difference in the initial model, w_lateral,n (τ

_j,l) is a trace for parameter n describing the weighting for each parameter sample where this weighting is a function of τ

_j,l which is a trace which at each parameter sample provides a measure of the local correlation between the traces j and l, #traces is the number of traces, #neighbors is the number of neighboring traces used in the calculation, #parameters is the number of parameters and L_p,lateral is the adjustable norm of said differences d_j,n,l.

48. The method according to claim 47, wherein the parameter difference d_j,l,n is defined as

(d_j,l,n)(t_k)=p_l,n(t_k+Δ

t_j,l,k)−

p_j,n(t_k)−

(p_initial,l,n(t_k+Δ

t_j,l,k)−

p_initial,j,n(t_k))

wherein Δ

t_j,l,k is the time shift at parameter sample k which time aligns the parameters of trace l to trace j at sample k, where surrounding trace samples are interpolated if at time t_k+Δ

t_j,k,l a sample is not defined, where r_j,l is now the local correlation incorporating the time shift and L_P,lateral is an adjustable norm on the parameter differences.

49. The method according to claim 37, wherein the seismic reflection data is determined from at least one of the following source-receiver combinations:

P-wave source and P-wave receiver, P-wave source and S-wave receiver, S-wave source and P-wave receiver, S-wave source and S-wave receiver.

50. The method according to claim 33, wherein a correction term is a measure for the differential time shifts between traces of measured reflection data stacks.

51. The method according to claim 50, wherein said measure comprises $F_{\mathrm{time}}=\sum _{}^{}$

i=2#⁢

 

⁢

stacks⁢

 

⁢

wtime,i⁢

∑

j=1#⁢

 

⁢

traces⁢

 

⁢

LP,time⁡

(τ

i,j-τ

0i,j)

wherein τ

_0,i,j is the trace with the initial time values of the time stretch and squeeze control points for stack i and trace j and τ

_i,j is the time of shifted control points, w_time,i is a weighting factor for stack i, #stacks is the number of stacks, #traces is the number of traces and L_P,time is an adjustable normalization factor of the difference between τ

_0,i,j and τ

_i,j.

52. The method according to claim 33, wherein a stabilization term is a measure for the parameter differences between reflection data acquisition surveys taken at different points in time.

53. The method according to claim 52, wherein said measure comprises $F_{\mathrm{timelapse}}=\sum _{}^{}$

k=2#⁢

 

⁢

surveys⁢

 

⁢

wsurvey,k⁢

∑

n=1#⁢

 

⁢

parameters⁢

 

⁢

wparameter,n⁢

∑

j=1#⁢

 

⁢

traces⁢

 

⁢

LP,timelapse⁡

(pj,n,k-pj,n,k-1)

wherein (p_j,n,k−

p<sub>j,n,k−

1</sub>) is the difference in parameters of trace j and parameter n between survey k and its time preceding survey k−

1, w_parameters,n is a weighting factor for each parameter n, w_survey,k is a weighting factor for each survey k, L_P,timelapse is an adjustable norm of said difference, #surveys represents the number of surveys and #parameters represents the number of parameters.

54. The method according to claim 33, wherein constraints and/or constraint functions are applied to one or more of the subsurface parameters;

and/or are applied to control changes of subsurface parameters between surveys taken on different times; and

/or constrain outside a specified subsurface zone the parameter changes between surveys to small values relative to changes expected within the specified subsurface zone; and

/or constrain the minimization by setting the subsurface parameters outside a specified subsurface zone from survey to survey at the same value.

55. The method according to claim 33, wherein for optimizing the objective function outside a specified subsurface zone only one or one set of different subsurface parameters are specified.

56. The method according to claim 33, comprising the generation of quality control information.

57. The method according to claim 56, wherein the quality information includes at least one of:

(i) synthetic data based on the optimized subsurface parameters;

(ii) residual data obtained by subtracting the synthetic data from the measured reflection data;

(iii) deviation data obtained by determining the deviations away from the initial subsurface parameters;

(iv) deviation data obtained by determining the deviations away from the corresponding functional relations; and

(v) deviation data obtained by determining the deviations away from well log data.

58. The method according to claim 33, wherein the reflection data is at least one of seismic data and time lapse data.

59. The method according to claim 58, wherein the time lapse data at each survey time comprises at least one seismic partial or full stack.

60. The method according to claim 33, wherein the subsurface parameters comprise elastic parameters.

61. The method according to claim 60, wherein the elastic parameters comprise pressure wave velocities and/or shear wave velocities and/or densities in the subsurface and/or the subsurface parameters comprise any mathematical relation between pressure wave velocities and/or shear wave velocities and/or densities.

62. The method according to claim 61, wherein the seismic data comprises at least two seismic partial or full stacks containing different angle dependant information on seismic reflections in the subsurface.

63. The method according to claim 33, wherein the subsurface parameters comprise compositional parameters representing the rock and fluid composition of the subsurface.

64. The method according to claim 33, wherein the reflection data is echo-acoustic data and the subsurface is human or mammal tissue or any other material.

65. A device for determining from measured reflection data on a plurality of trace positions one or more subsurface parameters, the device comprising:

(a) input means for inputting at least the measured reflection data and one or more initial subsurface parameters defining an initial subsurface model;

(b) processing means for;

(i) preprocessing the measured reflection data into a plurality of partial or full stacks;

(ii) specifying a wavelet or wavelet or wavelet field for each of the partial or full stacks of the measured reflection data;

(iii) calculating synthetic reflection data based on the specified wavelets or wavelet fields and the initial subsurface parameters; and

(iv) (a) optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data for a plurality of trace positions simultaneously;

or

(b) optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data for a plurality of trace positions simultaneously and one or more stabilization terms and/or one or more correction terms; and

(c) output means for outputting optimized one or more subsurface parameters.

66. A device for determining from measured reflection data on one or more trace positions a plurality of subsurface parameters, said device comprising:

(a) input means for inputting at least the measured reflection data and one or more initial subsurface parameters defining an initial subsurface model;

(b) processing means for;

(i) preprocessing the measured reflection data into a plurality of partial or full stacks;

(ii) specifying a wavelet or wavelet field for each of the partial or full stacks of the measured reflection data;

(iii) calculating synthetic reflection data based on the specified wavelets or wavelet fields and the initial subsurface parameters; and

(iv) (a) simultaneously optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data, for said plurality of subsurface parameters and for each trace position separately;

or

(b) simultaneously optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data for a plurality of subsurface parameters and one or more stabilization terms and/or one or more correction terms; and

(c) output means for outputting optimized subsurface parameters.

67. A method for determining from measured reflection data on a plurality of trace positions one or more subsurface parameters, said method comprising the steps of:

(a) preprocessing the measured reflection data into a plurality of partial or full stacks;

(b) specifying one or more initial subsurface parameters defining an initial subsurface model;

(c) specifying a wavelet or wavelet field for each of the partial or full stacks of the measured reflection data;

(d) calculating synthetic reflection data based on the specified wavelets and the initial subsurface parameters;

(e) optimizing an objective function, the objective function comprising the weighted difference between measured reflection data and synthetic reflection data, and one or more stabilization terms and one or more corrections terms, for a plurality of trace positions simultaneously; and

(f) outputting the optimized one or more subsurface parameters.

68. A device for determining from measured reflection data on a plurality of trace positions one or more subsurface parameters, said device comprising:

(a) input means for inputting at least the measured reflection data and one or more initial subsurface parameters defining an initial subsurface model;

(b) processing means for;

(i) preprocessing the measured reflection data into a plurality of partial or full stacks;

(ii) specifying a wavelet or wavelet field for each of the partial or full stacks of the measured reflection data;

(iii) calculating synthetic reflection data based on the specified wavelets or wavelet fields and the initial subsurface parameters; and

(iv) optimizing an objective function, comprising the weighted difference between measured reflection data and synthetic reflection data, and one or more stabilization terms and/or one or more correction terms for a plurality of trace positions simultaneously; and

(c) output means for outputting optimized one or more subsurface parameters.