GB2449497A

GB2449497A - Method and apparatus for processing electromagnetic response data

Info

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

Abstract

A method is provided for processing electromagnetic response data from a region of the earth. The method may be used for detecting subterranean conductivity anomalies, such as high-resistivity hydrocarbon reservoirs. Separate transverse electric response data and transverse magnetic response data are provided (15-23). Depth migration (24) is performed on the transverse electric and magnetic data using different depth migration extrapolation operators. The operators differ from each other in accordance with different properties of transverse electric and magnetic modes in an electromagnetically anisotropic medium.

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 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 GBO1/00419. Techniques of this type are also disclosed in Ellingsmd et al., 2002, "Remote sensing of hydrocarbon layers by seabed logging SBL: Results from a cmise offshore Angola", First Break, 21, No 10, 972-982, and in Eidesmo at elI 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 all, 2006, "Decomposition of electromagnetic fields into upgoing and downgoing components", Geophysics 71, No 5, G21 1-G223; Zhdanov et all, "Underground imaging by frequency-domain electromagnetic migration", Geophysics, 61, No 3, 666-682; Hokstad et all, 2006, "Anisotropic depth migration of marine controlled-source electromagnetic data", EAGE, 68th Ann. Intemat.

Mtg., Eur. Assoc. Geosc. Eng., Extended Abstracts A012; and Rosten et all, 2006, "3-D depth migration operators for marine controlled-source electromagnetic data", SEG, 76 Ann. Intemat. Mtg., Soc. Expl. Goephys., Expanded Abstract, 770-773.

According to a first aspect of the invention, there is provided a method of processing electromagnetic response data from a region of the earth for detecting subterranean conductivity or resistivity anomalies, comprising providing separate transverse electric response data and transverse magnetic response data, performing depth migration on the transverse electric data using a first depth migration extrapolation operator, and performing depth migration on the transverse magnetic data using a second depth migration extrapolation operator which differs from the first operator in accordance with different properties of transverse electric and magnetic modes in an electromagnetically anisotropic medium.

Such a technique is capable of detecting subterranean conductivity anomalies in regions of the earth, for example so as to detect, map or monitor subterranean hydrocarbon reservoirs. This technique makes no assumptions about electromagnetic isotropy and is therefore capable of being used to process data obtained from regions of the earth which are anisotropic in respect of their electromagnetic properties. The technique is responsive to regions having anisotropic conductivity but may also be used for monitoring or exploring regions having anisotropic permeability or anisotropic permittivity. Such techniques may therefore be used for electromagnetic data obtained from a wider range of geologies and to provide more reliable information about subterranean structure, such as the presence, shape and size of hydrocarbon reservoirs.

It has been found that it is possible to deal with anisotropy by processing transverse electric and magnetic response data using different operators. In particular, it is sufficient merely to select the appropriate different operators in order to take into account anisotropy of the medium which was explored or surveyed using electromagnetic techniques. Improved migration and imaging may therefore be performed irrespective of the presence and degree of electromagnetic anisotropy.

The transverse electric and magnetic data may be one-way response data. For example, the one-way response data may be upwardly-going response data.

The transverse electric and magnetic data may be quasi-data, that is, data that are only approximately TE or TM due to non-perfect preprocessing or anisotropy.

The transverse electric and magnetic data may be measured separately. As an alternative, the transverse electric and magnetic data may be derived from total field measurements.

The transverse electric and magnetic data may be frequency domain data. As an alternative, the method may comprise providing the transverse electric and magnetic data as time domain data and converting the time domain data to frequency domain data.

Each depth migration may comprise wavefield extrapolation and imaging. Each wavefield extrapolation may be performed in the frequency-wavenumber domain. The first and second operators may be phase shift operators.

The phase shift operators may be given by: 7&k. e

where y = i for down-going fields and y = -1 for migrated fields, and where: i = Pi is the imaginary unit; iz is the migration depth step; k for the first operator is given by 2(2 + k for the second operator is given by /K12 + y 2(1+ 211)(k + k and k are orthogonal horizontal wavenumbers; ic is the complex wavenumber in the ultralow frequency approximation 1iO4t01; * * i is given by 0 is the angular frequency; t0 is the permeability; and and a33 are the horizontal and vertical conductivities, respectively.

Each wavefield extrapolation may be performed in the frequency-space domain. The first and second operators may be discrete convolution filters. The convolution filters may be functions of the migration depth step, horizontal sampling interval, angular frequency and local conductivity. Coefficients of the convolution filters may be provided in the form of at least one pre-computed look-up table which is accessed during each wavefleld extrapolation.

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

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

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

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

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

According to a seventh aspect of the invention, there is provided output data produced by a method according to the first 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 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 I 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. 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. 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 interface 1 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 8, the angle of reflection (equal to the angle of incidence) by 9' and the transmitted angle by 92. The wavenumbers k1 and k1' in the medium 2 are equal to each other. The unit vector n is normal to the interface at the point of incidence. The medium 2 has electric permittivity Ci, magnetic permeability j.t' and conductivity o whereas the medium 3 has electric permittivity E2, magnetic permeability and conductivity cY2. (For simplicity, Figures 1 and 2 illustrate the isotropic case) Figure 3 illustrates an exploration vessel 5 at the surface of a water-column 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 11 is disposed on the sea bed 8. 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 two sets of steps, either 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, whereas the steps 20-23 are performed where the data are already separated into transverse electric and transverse magnetic components. A depth migration or extrapolation step 24 is common to both sets of steps.

In the step 15, the total field 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 measured transverse electric and transverse magnetic 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 up-going and down-going wavefields 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 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 24 may receive only one of the up- going and down-going wavefields. The step 24 performs depth migration of the up-going and/or down-going transverse electric and transverse magnetic components as described hereinafter. The output data from the step 24 may then be further processed or analysed in order to provide useful 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 backscattered from the air layer are contained in the TE mode.

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

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 TE and TM modes become different. The TE mode is sensitive only to the horizontal conductivity of the subsurface. The TM mode is 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 TE and TM migration operators are equal in the isotropic limit. In the proposed migration scheme, the depth migration operators for 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=i.tH, and VxH=JE, (1) Where E is the electric field, H is the magnetic field, p is the 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 if =cE, (2) where a is the conductivity tensor. Combining the equations above to eliminate the magnetic field gives the two-way diffusive wave equation JJE -a1E +ic4i0o11E1 =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, 1E1 = 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) jE =1oE =0, (4) since the divergence of a curl is identically zero. The conductivity tensor in a TI medium can be written as (5) where 0H and 033 are the horizontal and vertical conductivities, 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, aE, =-2fl(E +aE), (6) (7) 1 2q --where the electromagnetic anisotropy parameter is given by = O, (8) Substituting equation 6 for the divergence in equation 3 gives two coupled two-way equations for the horizontal components of the electric field (1+ 2q) + . E, + + 211,E. + i E, =0, (9) aE1. +(1+2fl)Er +aE), +icE, = 0, (10) where ic = .Jio4.t011 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 (1+21)(E: +,E:)+E: +icE: =0. (11) In the Fourier domain, equations 9 to 11 can be expressed as a 3 x 3 eigenvalue problem.

q2 -2qk 2flkrk,. 0 Er E, -2qk,k q2 -2rlk 0 E, =k E,. , (12) 0 0 q2 -2ri(k + k) E: E. where q2 = --k. The sign convention in the Fourier transform is such that -io and -* ik. Solving the characteristic equation gives three pairs of eigenvalues (13) where the positive and negative signs correspond to down-going and up-going plane waves, respectively. The three corresponding orthogonal and normalized eigenvectors can be written as = .-_ []anci = _L. [kr] and x = (14) where kr = jk + k,' is the radial horizontal wavenumber.

For a plane wave propagating with wavenumber k = (kr, k, k), from equation 14 it is clear that x is confined to the horizontal plane and that k = 0. Hence, the eigenvector corresponds to the TE mode with vertical wavenumber k?. The eigenvectors x2 and x belong to the degenerate eigenvalues k2 k3. Then, any linear combination y = a + is also an eigenvector with the same pair of eigenvalues. Since y x = o, x and x3 form an orthogonal basis for the TM mode with vertical wavenumber k? = k?. The derivation above has been given for downgoing and upgoing fields. The corresponding equations for the migrated field E;" 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 1-D background medium, standard one-way equations in the frequency-wavenumber (FK) domain for the down-going field E7 and migrated field E' are obtained as a..E,'=kE1', (15) where the vertical wavenumber k.. for the TE and TM modes is obtained from equation 13 as (16) where y = i for v = D and y = -1 for v = M. The solution to equation 15 is given by L(k,k,,,z+/z,O)=e"E(k,k,,,z,W), (17) which is the basis for the wavefield extrapolation step of FK migration.

The depth-stepping equation for E71 is numerically stable with exponential decay (like LID) and backward phase rotation (like E).

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

the electric field can be written as

L

E,'(x,v,z + Lz,o) = [Wv(,n&,n1w,ii,1)ET(x -tn&,y -niv,z,oj.

ni.rn-L In the following, equal spatial sampling distances are assumed in the two orthogonal horizontal directions, iy = Ex. The convolution operators W depend only on the normalized wavenumber i = ic4t, the anisotropy parameter i and the ratio z/x.

Hence, for a given L\z/Ax-ratio, the operator coefficients can be precompted for all relevant values of e and,, and stored in a look-up table. Computation of the finite impulse-response filter with complex-valued coefficients W'. (m&, nAv,i'1,i1) is posed as an inverse problem, minimizing the objective function in the L norm J=lI III" where W jtk, ,i is the discrete Fourier transform of W (m&, nEv, ,r) for a discrete set of horizontal wavenumbers iEk,1 and 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 W1D. The operator for W' 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 1 -D media, downward continuation is performed in the FK-domain by the phase-shift operator exp (ybzk..) 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 W' (in.x, n4v,i'1,11).

Figures 5a to 5c illustrate depth-stepping using these filter coefficients. The electric wavefield component E(m&,n4v,z) at angular frequency 0 is provided at each grid-point in a plane 30, as shown in Figure 5a. The filter coefficients W (m&, nv,1,11) are functions of the normalized (or dimensionless) wavenumbers.

The filter coefficients are pre-computed and made accessible in look-up table 31 where the respective depth-stepping length z, angular frequency and local conductivity is used to determine the correct operator at each grid point during the downward continuation process. In particular, as illustrated in Figure Sb, the filter look-up table index representing the filter coefficients W1(&,nv,i1,11) is accessed using the local grid-point conductivity, given the frequency and depth-stepping length as input.

Wavelength extrapolation of the local grid-point in the plane 30 from z to z + iz in the plane 32 by discrete convolution is performed according to equation 18 as illustrated in Figure 5c.

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 Ho [berg, 1988, mentioned hereinbefore.

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 K is defined by = Lx = (20) where / E [o, nI 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 [0,t 1' such that max K =nij=(l+i)7t. (21) Since /i = (i +1)! 1j the real and imaginary parts of K i are equal. Then the index I 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 /0toCii (22) icY 2 Introducing angular frequency o=2mf, (23) wheref is the linear frequency, the look-up table index I can be written as l=j7E, (24) where c is a precomputed function depending on the 3-D conductivity model, c(x,v, z)= ,1oQII(x,y,z). (25) In CSEM surveying, the maximum frequency is usually low. Typically, it is only necessary to consider the frequencies less than 10 Hz. Hence, / 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 factor s in the operator look-up table such that the scaled index L is given by L=s4J7c=sl, (26) and K1=LzKi (27) The maximum value of the scaled wavenumber corresponds to L = n. This gives = (1 * (28) Hence, in practice, the scaling implies that optimized convolution operators are computed on the interval [0,it Is].

The scaled index L given by equation 26 is valid for a fixed electromagnetic anisotropy parameter i. 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 t for different (discrete) values of i implying in general a 2D look-up table for wavefield extrapolation of TM and TE components in the FX-domain.

Claims

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

region of the earth for detecting subterranean conductivity or resistivity anomalies, comprising providing separate transverse electric response data and transverse magnetic response data, performing depth migration on the transverse electric data using a first depth migration extrapolation operator, and performing depth migration on the transverse magnetic data using a second depth migration extrapolation operator which differs from the first operator in accordance with different properties of transverse electric and magnetic modes in a electromagnetically anisotropic medium.
2. A method as claimed in claim 1, in which the transverse electric and magnetic data are one-way response data.
3. A method as claimed in claim 2, in which the one-way response data are upwardly-going response data.
4. A method as claimed in any one of the preceding claims, in which the transverse electric and magnetic data are quasi-data.
5. A method as claimed in any one of the preceding claims, in which the transverse electric and magnetic data are measured separately.
6. A method as claimed in any one of claims I to 4, in which the transverse electric and magnetic data are derived from total field measurements.
7. A method as claimed in any one of the preceding claimes, in which the transverse electric and magnetic data are frequency domain data.
8. A method as claimed in claim 7, comprising providing the transverse electric and magnetic data as time domain data and converting the time domain data to frequency domain data.
9. A method as claimed in any one of the preceding claims, in which each depth migration comprises wavefield extrapolation and imaging.
10. A method as claimed in claim 9, in which each wavefield extrapolation is performed in the frequency-wavenumber domain.
11. A method as claimed in claim 10, in which the first and second operators are phase shift operators.
12. A method as claimed in claim 11, in which the phase shift operators are given by: y&k.. e --

where y = i for down-going fields and y = -1 for migrated fields, and where: I = -/T is the imaginary unit; L\z is the migration depth step; k for the first operator is given by --y + k k for the second operator is given by.Iii2 +72(1 + 2)(k +k); k and k are orthogonal horizontal wavenumbers; i is the complex wavenumber.Jitxç.t0a11 * . C11 C33 i is given by 0.) is the angular frequency; is the permeability; and 0H and 033 are the horizontal and vertical conductivities, respectively.
13. A method as claimed in claim 9, in which each wavefield extrapolation is performed in the frequency-space domain.
14. A method as claimed in claim 13, in which the first and second operators are discrete convolution filters.
15. A method as claimed in claim 14, in which the convolution filters are functions of the migration depth step, horizontal sampling interval, angular frequency and local conductivity.
16. A method as claimed in claim 14 or 15, in which coefficients of the convolution filters are provided in the form of at least one precomputed look-up table, which is accessed during each wavefield extrapolation.
17. An apparatus arranged to perform a method as claimed in any one of the preceding claims.
18. A computer program for programming a computer to perform a method as claimed in any one of claims I to 16.
19. A computer containing a program as claimed in claim 18.
20. A computer-readable storage medium containing a program as claimed in claim 18.
21. Transmission across a network of a program as claimed in claim 18.
22. Output data produced by a method as claimed in any one of claims I to 16, by an apparatus as claimed in claim 17 or by a computer as claimed in claim 19.