GB2449509A

GB2449509A - Method for processing electromagnetic response data

Info

Publication number: GB2449509A
Application number: GB0717355A
Authority: GB
Inventors: Ketil Hokstad; Tage Roesten
Original assignee: Statoil ASA
Current assignee: Equinor ASA
Priority date: 2007-05-25
Filing date: 2007-09-07
Publication date: 2008-11-26
Also published as: GB2449497A; GB0717355D0; GB0710061D0

Abstract

A method for processing electromagnetic response data from a region of the earth maybe used for detecting subterranean conductivity anomalies, such as high-resistivity hydrocarbon reservoirs. A look-up table containing values of a frequency-space domain wavefield extrapolation operator addressed by horizontal spatial sampling intervals, extrapolation depth distance, a frequency parameter, and conductivity, is pre-computed and made available. For each of a plurality of points in the region of the earth, the look-up table is addressed by values of the horizontal spatial sampling intervals, extrapolation depth distance, the frequency parameter and local conductivity at the point. The appropriate operator value is retrieved from the look-up table and is used in performing wavefield extrapolation, for example as part of depth migration.

Description

Method of and an Apparatus for Processing Electromagnetic Response Data

The present invention relates to a method of and an apparatus for processing electromagnetic response data from a region of the earth, for example for detecting subterranean conductivity or resistivity anomalies. Such a method and apparatus may be used, for example, for processing land or marine electromagnetic surveying data for hydrocarbon exploration and for electromagnetic monitoring of hydrocarbon production.

A known technique for hydrocarbon exploration or hydrocarbon production monitoring below the sea bed is known as controlled source electromagnetic (CSEM) surveying or sea bed logging (SBL). Such a technique is disclosed for mapping hydrocarbons over shallow prospects in deep water in GBOI/00419. Techniques of this type are also disclosed in Ellingsrud et at., 2002, "Remote sensing of hydrocarbon layers by seabed logging SBL: Results from a cruise offshore Angola", First Break, 21, No 10, 972-982, and in Eidesmo et a!. 2002, "Sea Bed Logging (SBL), a new method for remote and direct identification of hydrocarbon filled layers in deepwater areas", The Leading Edge, 20, No 3, 144-152.

Various techniques are known for processing data gathered by SBL and other electromagnetic exploration techniques. Examples of processing techniques are disclosed in: GB 2411006; Amundsen et at., 2006, "Decomposition of electromagnetic fields into upgoing and downgoing components", Geophysics 71, No 5, G211-G223; Zhdanov et at., "Underground imaging by frequency-domain electromagnetic migration", Geophysics, 61, No 3, 666-682; Hokstad et al., 2006, "Anisotropic depth migration of marine controlled-source electromagnetic data", EAGE, 68th Ann. Intemat.

Mtg., Eur. Assoc. Geosc. Eng., Extended Abstracts, A012; R sten et al., 2006, "3-D depth migration operators for marine controlled-source electromagnetic data", SEG, 76th Ann. Intemat. Mtg., Soc. Expi. Goephys., Expanded Abstract, 770-773; and Hokstad et al., 2007, "3D anisotropic depth migration operators for marine controlled-source electromagnetic data", EAGE, 69th Ann. Intemat. Mtg., Eur. Assoc. Geosc. Eng., Extended Abstracts, D038.

According to a first aspect of the invention, there is provided a method of processing electromagnetic response data from a region of the earth, comprising the steps of: providing a look-up table containing values of an operator for wavefield extrapolation in the frequency-space domain addressed by horizontal spatial sampling intervals, extrapolation depth distance, a frequency parameter, and a local conductivity parameter; providing the electromagnetic response data; and, for each of a plurality of points in the region of the earth, addressing the look-up table by values of the horizontal spatial sampling intervals, the extrapolation depth distance, the frequency parameter, and the local conductivity parameter to retrieve one of the operator values and performing wavefield extrapolation using the retrieved operator value.

The look-up table may be further addressed by an electromagnetic anisotropy parameter.

The addressing step may comprise, for each of the points, addressing the look-up table by the values at the point.

The local conductivity parameter value used in the addressing step for each of at least some of the points may be the value at the point.

The local conductivity parameter value used in the addressing step for each of at least some of the points may be derived from the conductivity at at least one other point.

The providing step may comprise determining the operator values by minimising the objective function in the L norm: where wv (ik, jk. , ic, , q) is a discrete Fourier transform of W (mix, niy, K, , for a discrete set of horizontal wavenumbers iAk and jk, and ic and ii are a normalised wavenumber and an anisotropy parameter, respectively.

The response data may comprise frequency domain (x, y, z, co) data.

The response data may comprise electric and magnetic data.

The response data may comprise one-way response data.

The response data may comprise isotropic or transversely isotropic data.

The response data may comprise separate transverse electric response data and transverse magnetic response data. The transverse electric and magnetic data may be quasi-data.

The frequency parameter may be scaled frequency for increasing the sampling density of the look-up table at lower frequencies. The frequency parameter may comprise the product of a scaling factor and the square root of frequency.

The wavefield extrapolation may comprise depth extrapolation in the frequency-space domain.

The wavefield extrapolation may be followed by imaging to perform migration.

According to a second aspect of the invention, there is provided a method of preparing for processing of electromagnetic response data from a region of the earth, comprising: calculating values of a frequency- space domain wavefield extrapolation operator as a function of horizontal spatial sampling intervals, extrapolation depth distance, a frequency parameter, and a local conductivity parameter; ananging the values as a look-up table which is addressable by the horizontal spatial sampling intervals, the extrapolation depth distance, the frequency parameter, and the local conductivity parameter; and making the look-up table available for subsequent processing of response data.

The operator may be a discrete convolution filter.

According to a third aspect of the invention, there is provided an apparatus arranged to perform a method according to the first or second aspect of the invention.

According to a fourth aspect of the invention, there is provided a computer program for programming a computer to perform a method according to the first or second aspect of the invention.

According to a fifth aspect of the invention, there is provided a computer containing a program according to the fourth aspect of the invention.

According to a sixth aspect of the invention, there is provided a computer-readable storage medium containing a program according to the fourth aspect of the invention.

According to a seventh aspect of the invention, there is provided transmission across a network of a program according to the fourth aspect of the invention.

According to an eighth aspect of the invention, there is provided output data produced by a method according to the first or second aspect of the invention or by an apparatus according to the second aspect of the invention or by a computer according to the fourth aspect of the invention.

The invention will be further described, by way of example, with reference to the accompanying drawings, in which: Figures 1 and 2 are diagrams illustrating the propagation of transverse electric and

magnetic components, respectively, of a wavefield;

Figure 3 is a diagram illustrating an example of a sea bed logging technique for gathering electromagnetic data; Figure 4 is a flow diagram illustrating a method of processing electromagnetic response data constituting an embodiment of the invention; and Figures 5a to Sc are diagrams illustrating depth migration forming part of the method illustrated in Figure 4.

The method described hereinafter may be used to analyse any type of electromagnetic survey data acquired on land, in sea water or on the sea bed. For example, the data may be acquired by controlled-source techniques using active sources including a horizontal electric dipole source, a vertical electric dipole source, a horizontal magnetic dipole source, a vertical magnetic dipole source, a line transmitter, a circular transmitter, and any combination of these. The data may alternatively or additionally be acquired using a passive source such as the natural electromagnetic field of the earth (magnetotellurics). For marine surveying, the or each source need not be in the water-column or on the sea bed. For example, the or each source may be deployed anywhere in the subsurface, for instance in a borehole or a well or on land. The or each receiver may be in the water-column or on the sea bed but may also be deployed anywhere in the subsurface, such as in a borehole or a well or on land.

For the purpose of illustration only, the present technique will be described with reference to data obtained by SBL using a towed horizontal dipole source disposed in the water-column with the electromagnetic receivers disposed as an array on the sea bed.

The present technique provides at least partially separate processing of the electric and magnetic data into one-way response data. The data may be isotropic or transversely isotropic. For example, the present technique may provide at least partially separate processing of the transverse electric mode response data and the transverse magnetic mode response data. Figure 1 illustrates the propagation of the transverse electric mode at an interface 1 between a first medium 2 and a second medium 3 of different electromagnetic properties. For the transverse electric mode, the magnetic field vector H is in the plane of incidence and the electric field vector E is in the horizontal plane perpendicular to the plane of incidence. The transverse magnetic mode is illustrated in Figure 2 with the electric field vector E in the plane of incidence and the magnetic field vector H in the horizontal plane perpendicular to the plane of incidence. The wavenumber vector k, the electric field vector E and the magnetic field vector H are shown for a right-handed coordinate system. The incoming fields incident at the interfice I have the subscript "1", the reflected fields have the subscript "1", and the transmitted fields have the subscript "2". The angle of incidence is denoted by 9, the angle of reflection (equal to the angle of incidence) by 0' and the transmitted angle by 02. The wavenumber is k1 in the medium I and k2 in the medium 2. The unit vector n is normal to the interface at the point of incidence. The medium 2 has electric permittivity ci, magnetic permeability iti and conductivity a whereas the medium 3 has electric pennittivity C2, magnetic permeability lt2 and conductivity 2. (For simplicity, Figures 1 and 2 illustrate the isotropic case) Figure 3 illustrates an exploration vessel 5 at the surface of a water-colunm comprising a region of sea 6. The water-column 6 lies above a subterranean region 7 of the earth with the sea bed 8 forming an interface between the sea and the earth. A subterranean feature is illustrated at 9 and comprises a region or formation of relatively high electrical resistivity. Such a region may comprise a hydrocarbon reservoir which is to be detected, mapped or monitored. However, the present technique may be used to explore other structures of high resistivity (low electrical conductivity), such as salt or basalt formations so as to obtain information about the geological structure of the region of the earth of interest.

By way of example, the vessel 5 tows a horizontal electric dipole source 10, which operates as a controlled source in the known SBL technique. A two-dimensional array of receivers or detectors Ills disposed on the sea bed 8. Note however that the receiver layout can be arbitrary. As the source 10 is towed, it is continuously or repeatedly actuated and the measurements performed by the array of detectors 11 are returned to the vessel 5 and recorded. The recorded data are subsequently processed by the present technique.

Figure 4 illustrates three sets of steps, any one of which may be performed depending on the nature of the data obtained by the SBL exploration. Thus, the steps 15-19 are performed where the electromagnetic data represent total field data without any separation of modes, the steps 20-23 are performed where the data are already separated into up-and down-going wavefield components, and the steps 24-27 are performed where the data are already separated into transverse electric and transverse magnetic components. A depth migration or extrapolation step 28 is common to all three sets of steps.

In the step 15, the total afield electric and magnetic field data are measured during the SBL survey and made available to a suitable processing system, such as a programmed computer or array of computers. In a step 16, the data are rotated so as to be in a desired Cartesian coordinate system. Techniques for performing such rotation are well-known are will not be described further.

In the step 17, the data are transformed from the time domain to the frequency domain, for example by a numerical Fourier transform. Again, such transforms are well-known and will not be described further. In the step 18, the data are processed so as to separate the up-going and down-going wavefields. Such separation techniques are also known and an example is disclosed in the publication by Amundsen et al, 2006, mentioned hereinbefore. In the step 19, the data are separated into transverse electric (TE) and transverse magnetic (TM) modes. Techniques for performing the separation are known and an example is disclosed in GB24 11006.

The order of the steps 16 to 19 is given by way of example. These steps may be performed in any desired or convenient order and some steps are optional.

The step 20 illustrates the case where the SBL exploration provides separate up-and down-going electric and magnetic wavefield components. The step 21 performs coordinate rotation and may use the same techniques as the step 16. The step 22 performs the transform from the time domain to the frequency domain and may use the same techniques as the step 17. The step 23 separates the transverse electric and transverse magnetic components and may use the same techniques as the step 19.

Again, the steps 21 to 23 may be performed in any desired or convenient order and some steps are optional.

The step 24 illustrates the case where the SBL exploration provides separate measured transverse electric and transverse magnetic components. The step 25 performs coordinate rotation and may use the same techniques as the step 16. The step 26 performs the transform from the time domain to the frequency domain and may use the same techniques as the step 17. The step 27 separates the up-going and down-going wavefields and may use the same techniques as the step 18. Again, the steps 25 to 27 may be performed in any desired or convenient order and some steps are optional.

The step 28 receives the separated up-going and down-going wavefields in the frequency domain with respect to the desired coordinate system and with the transverse electric and transverse magnetic components separated. Alternatively, the step 28 may receive only one of the up- going and down-going wavefields. The step 28 performs depth migration of the up-going and/or down-going wavefields or the up-going and/or down-going transverse electric and transverse magnetic components as described hereinafter. The output data from the step 28 may then be further processed or analysed in order to provide useflil information about the nature of the structure 9.

The TE and TM modes have significantly different geophysical properties: 1. The TE mode is mostly sensitive to near-surface and overburden properties of the subsurface. Also, the electric and magnetic fields backscauered from the air layer are mostly contained in the TE mode.

2. The TM mode is mostly sensitive to thin high resistive objects in the subsurface, such as hydrocarbon-filled reservoirs.

Also, the TE and TM modes react differently to anisotropy, which is frequently observed in a stratified subsurface.

For the purpose of hydrocarbon exploration, it is therefore advantageous to perform separation of the TE and TM modes at an early stage of the data processing.

Subsequently, appropriate processing can be applied separately to the TE and TM modes, to extract the specific information contained in the respective modes.

If the conductivity (or its inverse, the resistivity) of the subsurface is anisotropic, the propagation and scattering properties of the electric and magnetic data and including the TE and TM modes become different. The isotropic electric and magnetic data and the TE mode ar sensitive only to the horizontal conductivity of the subsurface. The anisotropic electric and magnetic data and the TM mode are sensitive to both horizontal and vertical conductivity. 3D FK and FD migration operators are described for depth extrapolation of TE and TM modes in transversely isotropic media. The electric and magnetic migration operators and the TE and TM migration operators are equal in the isotropic limit. In the proposed migration scheme, the depth migration operators for electric and magnetic data and for the TE and TM modes can be used in combination with a non-local imaging condition, accounting for guided waves. In general anisotropic media (beyond vertical transverse isotropic (TI)), the decomposition into TE and TM modes is not exact and the terms quasi-TM and quasi-TE are used to describe this case.

In the ultra-low frequency approximation, neglecting displacement currents, the Faraday and Ampere equations in the frequency domain can be approximated by VXE=jO!LH, and VXH=JE, (1) Where E is the electric field, H is the magnetic field, I = i:i is the imaginary unit, t is the magnetic permeability, and o is the angular frequency. Without external sources, the current density jE in a linear anisotropic medium is given by Ohm's law (using Einstein's sum notation) jE =oE (2) where o is a component of the conductivity tensor. Combining the equations above to eliminate the magnetic field gives the two-way diffusive wave equation aJaJE, -a1a 1E1 +io.uT,.E1 = 0. (3) The medium is assumed to be smooth such that the spatial derivatives of the conductivity can be neglected. This is justified by the fact that the goal is to perform migration, where the electric field is propagated in a smooth background medium.

In the isotropic case, aE1 = 0 and equation 3 reduces to the scalar diffusive Helmholz equation for each component separately. This does not hold in the anisotropic case.

However, from Ampere's law (with zero external source current) ajE =ajaq,Ej =0, (4) since the divergence of a curl is identically zero. The conductivity tensor in a TI vertical (TIV) medium can be written as (Ti' 0= a,, , (5) where 0jj and a are the horizontal and vertical conduct ivities, respectively.

Substituting the TI conductivity tensor above in equation 4, the divergence of the electric field can be expressed in two alternative ways which will both be used below, I0 a,E, =-2q(E +E), (6) = (7) where the electromagnetic anisotropy parameter is given by 1011-033, (8) Substituting equation 6 for the divergence in equation 3 gives two coupled two-way equations for the horizontal components of the electric field (l 211)aE +E aE 2TlaaE1 +1cE =0, (9) +(I+2r)E +aE. 2flaxa1E, +1c12E =0, (10) where ic, = ,Jioq.ta,, is the complex wavenumber. Returning to equation 3 and using equation 7 to eliminate the divergence of the electric field gives a separate two-way equation for the vertical component (l+2q)(2E +aE..)+a2E +icE.. =0. (11) In the Fourier domain, equations 9 to 11 can be expressed as a 3 x 3 eigenvalue problem.

q2 -2r1k -2flkXk). 0 E -2T)kxk). q2 -2rk 0 E. =k E, (12) 0 0 -2i (k + k) E E. Where k and k,, are the horizontal spatial wavenumbers, and q2 = ic' -k -k. The sign convention in the Fourier transform is such that 3, <--1w and -1k1. Solving the characteristic equation gives three pairs of eigenvalues (vertical wavenumbers) k -(k + k), and k2 = k3 = /ic -(1 + 2r)(k + k), (13) where the positive and negative signs correspond to down-going and up-going plane waves, respectively. The three corresponding orthogonal and normalized eigenvectors canbewrittenas -k. k 0 1c and x(2)=_L Ic,, and 0, (14) 0 ro 1 where Ic, = Jk + k is the radial horizontal wavenumber.

For a plane wave propagating with wavenumber Ic = (, k, k), from equation 14 it is clear that is confined to the horizontal plane and that k -x(1) = 0. Hence, the eigenvector xW corresponds to the TE mode with vertical wavenumber k!'. The eigenvectors x2 and belong to the degenerated eigenvalues k2 = Then, any linear combination y = a x'2 + f3x3' is also an eigenvector with the same pair of eigenvalues. Since y =0, and x form an orthogonal basis for the TM mode with vertical wavenumber /d2 = k!3. The derivation above has been given for down-going and up-going fields. The corresponding equations for the migrated field E,M are derived in the same way, using the adjoint two-way equation which is obtained by changing the sign of the last term in equations 9 to 11.

In a I -D background -medium, standard one-way equations in the frequency-wavenumber (FK) domain for the down-going field E and migrated field E' are obtained as (15) where the vertical wavenumber /c for the TE and TM modes is obtained from equation 13 as k!'E +y2(k+k2jand k!TM =.I1c+y2(I+2q)(k2+k2) (16) where = i for v = D (down-going field) and y = -1 for v M (migrated field) The solution to equation 15 is given by +&,co) = eTE'(k,k.,z,W), (17) which is the basis for the wavefleld extrapolation step of FK migration.

The depth-stepping equation for E7' (migrated field) is numerically stable with exponential decay (like E; down-going field) and backward phase rotation (like E;

up-going field).

The FK migration operators are accurate up to 90 degrees from the vertical, but limited to 1-D background media. To relax the l-D background assumption, the phase-shift operator in equation 17 may be replaced by discrete convolution filters in. the frequency-space (FX) domain, as proposed for seismic data by Holberg, 1988, "Towards optimal one-way wave propagation", Geophsical Prospecting, 36, 99-114.

The extrapolation of the electric field can be written as

L

Ej"(x,y,z + &,o) = [Wt(mAx,nty,ie, ,rl)E,t (x -m&,y -n4v,z,co]. (18) m.n--L In the following, without loss of generality, equal spatial sampling distances are assumed in the two orthogonal horizontal directions, iy ix. The convolution operators W" depend only on the normalized wavenumber ie = 1cx, the anisotropy parameter ij and the ratio AzJAx. Hence, for a given zJx-ratio, the operator coefficients can be precompted for all relevant values of ie and i and stored in a look-up table. Computation of the finite impulse-response filter with complex-valued coefficients W'(m,nv,1,r) is posed as an inverse problem, minimizing the objective function in the L norm J=II W'(iMc1,r) -e7 where wv (ii.k, jik,, , i,r) is the discrete Fourier transform of W (mtSx, n4y, ie,q) for a discrete set of horizontal wavenumbers itilç and f/i/c,. The dispersion relation for diffusive EM fields is smooth and continuous for all wavenumbers. Hence, it is not necessary to introduce a dip-limitation on the corresponding 3-D filter operators to obtain a stable depth migration scheme for CSEM data. The optimization is generally performed for all wavenumbers and the real and imaginary parts of the filter operator are optimized separately. In practice, it is necessary to compute tabulated filter coefficients only for W''. The operator for WM may then be obtained by complex conjugation.

Depth migration using one-way equations comprises two steps:

Wavefield extrapolation

* Imaging Wavefield extrapolation used in depth-stepping migration is sometimes called downward continuation. In l-D media, downward continuation is performed in the FK-domain by the phase-shift operator exp (y/izk..) as described in equation 17. In 3-D media, on the other hand, downward continuation can be performed by space-variant convolution in the FX-domain as described in equation 18 using the optimized filter coefficients WV (,n/ix, iz4v,ie1,q).

Figures 5a to 5c illustrate depth-stepping using these filter coefficients. The electric wavefield component E' (rn/ix, ntiy, z,0) at angular frequency o is provided at each grid-point in a plane 30, as shown in Figure 5a. The filter coefficients W"(nthx,nAv,i1,1) are functions of the normalized (or dimensionless) wavenumber and anisotropy parameter. The filter coefficients are pre-computed and made accessible in look-up table 31 where the respective horizontal spatial sampling intervals, depth-stepping length, angular frequency and local conductivity (x, y, z) is used to determine the correct operator at each grid point during the downward continuation process. In particular, as illustrated in Figure 5b, the filter look-up table index representing the filter coefficients Wv(rnIXr,náy,x1,11) is accessed using the local grid-point conductivity, given the frequency, horizontal spatial sampling intervals, depthstepping length, and anisotropy parameter as input. Wavelength extrapolation of the local grid-point in the plane 30 from z to z + tz in the plane 32 by discrete convolution is performed according to equation 18 as illustrated in Figure Sc.

Once the filter operators are designed, the resulting wave propagation algorithm is simple and is eminently suited for implementation on large parallel computing systems, for example as disclosed in the publication by Holberg, 1988, mentioned hezeinbefore.

The FK and FD depth-migration methods based on one-way equations do not handle wavefield amplitudes accurately. In fact, these schemes were never assumed to do so. Hence, the result of depth migration is a relative structural image of the resistivity contrasts. To compute estimates of the subsurface parameter model, full inversion must be used.

The normalized wavenumber Wi is defined by 25. = Jioit0c,& = !, (20) where / E [o, n} is the index in the operator look-up table and n is the maximum index in the table corresponding to the normalized Nyquist wavenumber. The real and imaginary parts of the normalized wavenumber are defined on the interval [o,], such that max x1 =n1Ci=(1+j)it. (21) Since.11= (1+ i)i the real and imaginary parts of ici are equal. Then the index F can be computed using only the real part (or only the imaginary part) of the normalized wavenumber. Comparing equations 20 and 21, the index I can be expressed as i=!!4 I0ui. (22) itY 2 Introducing angular frequency = 2mf, (23) wheref is the linear frequency, the look-up table index (can be written as l=Jje, (24) where c is a precomputed function depending on the 3-D conductivity model, (25) In CSEM surveying, the maximum frequency is usually low. Hence, I will always be an index in the lower part of the look-up table. On the other hand, for improved accuracy it is desirable to have a dense sampling of the part of the look-up table that is frequently used. A convenient way to achieve this is to introduce a scaling factors in the operator look-up table such that the scaled index L is given by L = = (26) and Ki=LiK (27) The maximum value of the scaled wavenumber corresponds to L = n. This gives malc K1 =n!K1=t (28)

S

Hence, in practice, the scaling implies that optimized convolution operators are computed on the interval [o,, Is].

The scaled index L given by equation 26 is valid for a fixed electromagnetic anisotropy parameter q. For the isotropic case, the look-up table will be 1-D. In the case of anisotropic media, it is, in general, necessary to calculate filter operators for (discrete) values of L between 0 and r for different (discrete) values of Ti implying in general a 2-D look-up table for wavefield extrapolation of isotropic and anisotropic electric and magnetic data and TM and TE components in the FX-domain.

Claims

CLAIMS: 1. A method of processing electromagnetic response data from a

region of the earth, comprising the steps of: providing a look-up table containing values of an operator for wavefield extrapolation in the frequency-space domain addressed by horizontal spatial sampling intervals, extrapolation depth distance, a frequency parameter, and a local conductivity parameter; providing the electromagnetic response data; and, for each of a plurality of points in the region of the earth, addressing the look-up table by values of the horizontal spatial sampling intervals, the extrapolation depth distance, the frequency parameter, and the local conductivity parameter to retrieve one of the operator values and performing wavefield extrapolation using the retrieved operator value.
2. A method as claimed in claim 1, in which the look-up table is further addressed by an electromagnetic anisotropy parameter.
3. A method as claimed in claim I or 2, in which the addressing step comprises, for each of the points, addressing the look-up table by the values at the point.
4. A method as claimed in any one of the preceding claims, in which the local conductivity parameter value used in the addressing step for each of at least some of the points is the value at the point.
5. A method as claimed in any one of the preceding claims, in which the local conductivity parameter value used in the addressing step for each of at least some of the points is derived from the conductivity at at least one other point.
6. A method as claimed in any one of the preceding claims, in which the providing step comprises determining the operator values by minimising the objective function in the L norm: =I1w(kx,ity,, ,n)-eII" where W"(ik,q) is a discrete Fourier transform of W (mix, niy, K, , q) for a discrete set of horizontal wavenumbers iL.k and jik, and ic and i are a normalised wavenumber and an anisotropy parameter, respectively.
7. A method as claimed in any one of the preceding claims, in which the response data comprise frequency domain data.
8. A method as claimed in any one of the preceding claims, in which the response data comprise electric and magnetic data.
9. A method as claimed in any one of the preceding claims, in which the response data comprise one-way response data.
10. A method as claimed in any one of the preceding claims, in which the response data comprise isotropic or anisotropic data.
11. A method as claimed in any one of the preceding claims, in which the response data comprise separate transverse electric response data and transverse magnetic response data.
12. A method as claimed in claim 11, in which the transverse electric and magnetic data are quasi-data.
13. A method as claimed in any one of the preceding claims, in which the frequency parameter is scaled frequency for increasing the sampling density of the look-up table at lower frequencies.
14. A method as claimed in claim 13, in which the frequency parameter comprises the product of a scaling factor and the square root of frequency.
15. A method as claimed in any one of the preceding claims, in which the wavefield extrapolation comprises depth extrapolation in the frequency-space domain.
16. A method as claimed in any one of the preceding claims, in which the wavefield extrapolation is followed by imaging to perform migration.
17. A method of preparing for processing of electromagnetic response data from a region of the earth, comprising: calculating values of a frequency-space domain wavefield extrapolation operator as a function of horizontal spatial sampling intervals, extrapolation depth distance, a frequency parameter, and a local conductivity parameter; arranging the values as a look-up table which is addressable by the horizontal spatial sampling intervals, the extrapolation depth distance, the frequency parameter, and the local conductivity parameter; and making the look-up table available for subsequent processing of response data.
18. A method as claimed in any one of the preceding claims, in which the operator is a discrete convolution filter.
19. An apparatus arranged to perform a method as claimed in any one of the preceding claims.
20. A computer program for programming a computer to perform a method as claimed in any one of claims 1 to 18.
21. A computer containing a program as claimed in claim 20.
22. A computer-readable storage medium containing a program as claimed in claim 20.
23. Transmission across a network of a program as claimed in claim 20.
24. Output data produced by a method as claimed in any one of claims I to 18, by an apparatus as claimed in claim 19 or by a computer as claimed in claim 21.