Modeling gravity and tensor gravity data using poisson's equation for airborne, surface and borehole applications

CROSS REFERENCE TO RELATED APPLICATIONS

0001This application is a Continuation-in-Part of U.S. patent App. Ser. No. 10/236,204 filed on Sep. 6, 2002 now U.S. Pat. No.6,675,097. App. Ser. No. 10/236,204 claims priority from U.S. Provisional App. Ser. No. 60/318,083 filed on Sep. 7, 2001. App. Ser. No. 10/236,204 is a Continuation-in-Part of U.S. patent app. Ser. No. 09/580,863 filed on May 30, 2000, now U.S Pat. 6,502,037. App. Ser. No. 09/580,863 is a Continuation-in-Part of U.S. patent app. Ser. No. 09/285,570 filed on Apr. 2, 1999, now U.S. Pat. No. 6,278,948; app. Ser. No. 09/399,218 filed on Sep. 17, 1999, now U.S. Pat. No. 6,424,918; and U.S. patent app. Ser. No. 09/405,850 filed on Sep. 24, 1999, now U.S. Pat. No. 6,430,507. App. Ser. No. 09/399,218 filed on Sep. 17, 1999, now U.S. Pat. No. 6,424,918 and U.S. patent app. Ser. No. 09/405,850 filed on Sep. 24, 1999, now U.S. Pat. No. 6,430,507 are both Continuation-in-Part applications of U.S. patent app. Ser. No. 09/285,570 filed on Apr. 2, 1999, now U.S. Pat. No. 6,278,948.

BACKGROUND OF THE INVENTION

00021. Field of the Invention

0003The present invention pertains to processing potential fields data using vector and tensor data along with seismic data and more particularly to the modeling and inversion of gravity data which may also be combined with seismic data for hydrocarbon exploration and development.

00042. Related Prior Art

0005Exploration for hydrocarbons in subsurface environments containing anomalous density variations has always presented problems for traditional seismic imaging techniques by concealing geologic structures beneath zones of anomalous density. Many methods for determining parameters related subsurface structures and the extent of the highly anomalous density zones exist.

0006U.S. Pat. No. 4,987,561, titled “Seismic Imaging of Steeply Dipping Geologic Interfaces, issued to David W. Bell, provides an excellent method for determining the side boundary of a highly anomalous density zone. This patent locates and identifies steeply dipping subsurfaces from seismic reflection data by first identifying select data which have characteristics indicating that seismic pulses have been reflected from both a substantially horizontal and a steeply dipping interface. These data are analyzed and processed to locate the steeply dipping interface. The processed data are displayed to illustrate the location and dip of the interface. This patent, while helping locate the boundaries, provides nothing to identify the subsurface formations on both sides of the boundary.

0007There have also been methods for identifying subsurface formations beneath anomalous zones using only seismic data to create a model and processing the data to identify formations in light of the model. By further processing reflection seismic data, the original model is modified or adjusted to more closely approximate reality.

0008An example of further processing seismic data to improve a model is U.S. Pat. No. 4,964,103, titled “Three Dimensional Before Stack Depth Migration of Two Dimensional or Three Dimensional Data,” issued to James H. Johnson. This patent provides a method of creating a three-dimensional model from two dimensional seismic data. This is done by providing a method of ray tracing to determine where to move trace segments prior to stacking. The trace segments are scaled to depth, binned, stacked and compared to the seismic model. The model can then be changed to match the depth trace segments that will be stacked better, moved closer to their correct three-dimensional position and will compare better to the model. This patent uses a computationally extensive seismic process to modify a seismic model that may be inaccurate.

0009One source of geologic exploration data that has not been used extensively in the past is potential fields data, such as gravity and magnetic data, both vector and tensor data and using potential fields data in combination with seismic data to provide a more accurate depth model or to derive a velocity model.

0010Gravity gradiometry has been in existence for many years although the more sophisticated versions have been held as military secrets. The measurement of gravity became more usable in the late eighteen hundreds as measuring instruments with greater sensitivity were developed. Prior to this time, while gravity could be measured, the gravity gradient due to nearby objects could not be reliably measured.

0011It has been known since the time of Sir Isaac Newton that bodies having mass exert a force on each other. The measurement of this force can identify large objects having a change in density even though the object is buried beneath the earth's surface or in other ways out of sight.

0012Exploration for hydrocarbons in subsurface environments containing anomalous density variations such as salt formations, shale diapers and high pressure zones create havoc on seismic imaging techniques by concealing geologic structures beneath zones that are difficult to image. Often the imaging problem is due to velocity anomalies directly related to the anomalous densities of these structures. By utilizing gravity, magnetic and tensor gravity field measurements along with a robust inversion process, these anomalous density zones can be modeled. The spatial resolution obtained from this process is normally lower than that obtained from reflection seismic data. However, models obtained from gravity and magnetic data can provide a more accurate models for the seismic processing. Using the potential fields data models as an aid in seismic depth imaging processing greatly enhances the probability of mapping these concealed geologic structures beneath the zones of anomalous density.

0013Zones of anomalous density may also be associated with zones of anomalous fluid pressure. Typically, while drilling an oil or gas well, the density of the drilling mud must be controlled so that its hydrostatic pressure is not less than the pore fluid pressure in any formation along the uncased borehole. Otherwise, formation fluid may flow into the wellbore, and cause a “kick.” Kicks can lead to blowouts if the flow is not stopped before the formation fluid reaches the top of the well. If the fluid contains hydrocarbons, there is a serious risk of an explosion triggered by a spark. For this reason, wellbores are drilled with a slight excess of the borehole fluid pressure over the formation fluid pressure.

0014A large excess of the borehole fluid pressure over the formation fluid pressure, on the other hand, is also undesirable. Fractures in the borehole wall may result in loss of circulation of the drilling fluid, resulting in stuck drill strings, time delays and greater costs. Serious formation damage may also occur that can decrease the amount of recoverable minerals.

0015Pressure prediction is done by estimating certain key parameters that include the overburden stress or confining stress, which is defined as the total lithostatic load on a rock volume at a given depth, and the effective stress, which is defined as the net load on the grain framework of the rock at a given depth. These two relations are then used in the Terzaghi effective stress law to estimate the fluid or pore pressure. Terzaghi's law states that: <br /><i>Pc=Pe+Pp </i>where:<ul id="ul0001" list-style="none"><li id="ul0001-0001" num="0000"><ul id="ul0002" list-style="none"><li id="ul0002-0001" num="0016">(Pc)=the confining stress</li><li id="ul0002-0002" num="0017">(Pe)=the stresses born by the grains, and</li><li id="ul0002-0003" num="0018">(Pp)=the stress born by the fluid.</li></ul></li></ul>

0019Some workers treat a special case of Terzaghi's law where the confining stress is assumed to be the mean stress as opposed to the vertical confining stress. It should be acknowledged that this difference exists, but that it does not effect the embodiments of the present invention as they will pertain to estimating the total overburden load, which can then be converted to either vertical confining stress or mean stress based on the stress state assumptions that are made. A prior art method for estimating confining stress is to use a density log from a nearby calibration well and integrate the density data to obtain the overburden load. This calibration is applied from the mudline down to depths usually beyond the depth of sampling to predict the overburden away from the calibration well.

0020It has long been recognized that velocities of seismic waves through sedimentary formations are a function of “effective stress,” defined as the difference between the stress caused by the overburden and the pore fluid pressure. A number of methods have been used to measure the seismic velocities through underground formations and make an estimate of the formation fluid pressure from the measured velocities. Plumley (1980) and U.S. Pat. No. 5,200,929 issued to Bowers, (the '929 patent) describe a method for estimating the pore fluid pressure at a specified location. The method also accounts for possible hysteresis effects due to unloading of the rock formation. The method utilized a pair of sonic velocity-effective stress relations. One relationship is for formations in which the current effective stress is the highest ever experienced. A second relationship is used when the effective stress has been reduced from the maximum effective stress experienced by the rock and hysteresis must be accounted for.

0021The '929 patent uses density data from nearby wells or from a geologically similar well to obtain the overburden stress. In most circumstances, the overburden stress may be adequately described by general compaction models in which the density increases with depth, giving rise to a corresponding relation for the relation between depth and overburden. In the absence of well control, determination of the overburden stress even within a sedimentary column is problematic. Furthermore, there are circumstances in which the model of a density that increases uniformly with depth is not valid. In such cases, the assumption of increasing density with depth is violated and a different approach to estimation of the overburden stress is needed.

0022There are several types of situations that may arise wherein a model of density increasing with depth and compaction is not valid. In the first case, there is a region of abnormally high density in the subsurface, usually of magmatic origin. The region could consist of an extrusive or intrusive volcanic material having relative density of 2.8 or higher. When such a formation is present within a sedimentary section where the relative density is typically between 2.4 and 2.65, the result is an increase in the overburden stress underneath the formation over what would be determined by prior art calculations. On the other hand, a region of abnormally low density may occur from salt bodies (2.10) or shale diapirs. In such a case, the overburden stress is abnormally low compared to what would be determined by prior art methods. In either case, even if the effective stress could be determined from seismic velocity measurements, a formation fluid pressure determination based on a prior art density model would be invalid.

0023<figref idref="DRAWINGS">FIG. 1</figref> is a seismic section illustrating an area of interest having highly anomalous density zones such as salt domes. A rounded interface can be seen as masking the formations below. For this data set, the lower boundary cannot be determined by normal seismic data processing methods.

0024There is a need for a method to image subterranean formations which are responsive to potential fields and non-potential fields data. There is a need for a method to combine these data types to extract more useful information than either data type has provided before. There is a need for a method that can more accurately determine the density of the subsurface in three dimensions away from and deeper than the limits of density from well control. The present invention satisfies this need.

SUMMARY OF THE INVENTION

0025The present invention is a method for determining a parameter of interest of a region of interest of the earth. At least one component of potential fields data is measured at a plurality of locations over a region of interest including a subterranean formation of interest. The potential field data are selected from magnetic data and gravity data. An initial geophysical model is determined for the region including the subterranean formation of interest. For the model, geophysical tensor data is updated using a forward model at a plurality of locations using a High Order Compact Finite Difference method and may include conjugate gradient methods. The method may be used to determine potential field measurements within a source region that is difficult to formulate using standard intregal equation approach where the Green's function is evaluated. Green's functions evaluations are expensive and require significant memory to store, where as in this formulation only requires storing a sparse stencil. A difference between the estimated model value and the measured value of the potential field measurements are determined at the plurality of locations. Based on this difference the geophysical model is updated. The model is iteratively updated and compared to the measured data until the differences reach an acceptable level and the parameter of interest has been determined.

BRIEF DESCRIPTION OF THE DRAWINGS

0026The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawings(s) will be provided by the Office upon request and payment of the necessary fee. The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

0027<figref idref="DRAWINGS">FIG. 1</figref> is a seismic section of an area having anomalous density zones such as a salt dome;

0028<figref idref="DRAWINGS">FIG. 2A</figref> illustrates a 19 point finite difference stencil with matrix coefficients;

0029<figref idref="DRAWINGS">FIG. 2B</figref> illustrates the finite difference stencil of <figref idref="DRAWINGS">FIG. 2A</figref> displayed as a cube;

0030<figref idref="DRAWINGS">FIG. 3</figref> illustrates a comparison of an analytic solution with numerically calculated model data;

0031<figref idref="DRAWINGS">FIG. 4</figref> illustrates a comparison of an integral formulation to a differential formulation;

0032<figref idref="DRAWINGS">FIG. 5A</figref> illustrates a realistic model with parameters which approximate a salt diapir;

0033<figref idref="DRAWINGS">FIG. 5B</figref> illustrates the finite difference forward modeling results of a realistic model with parameters which approximate a salt diapir of <figref idref="DRAWINGS">FIG. 5A</figref>;

0034<figref idref="DRAWINGS">FIG. 6</figref> illustrates borehole gravity apparent density depth dependence;

0035<figref idref="DRAWINGS">FIG. 7</figref> illustrates that the calculation of apparent density from borehole gravity using only the vertical component of the field results in an anomalously high density value;

0036<figref idref="DRAWINGS">FIG. 8</figref> illustrates the contribution of each tensor component in the apparent density in borehole measurements;

0037<figref idref="DRAWINGS">FIG. 9</figref> is a flow chart illustrating geophysical model development;

0038<figref idref="DRAWINGS">FIG. 10</figref> is flow chart of an embodiment of the invention using constrained optimization;

0039<figref idref="DRAWINGS">FIG. 11</figref> is flow chart of an embodiment of the invention using Laplace's Equation and/or the Equivalent Source method; and

0040<figref idref="DRAWINGS">FIG. 12</figref> is flow chart of an embodiment of the invention using a Poisson's Equation method solution.

DESCRIPTION OF THE PREFERRED EMBODIMENT

0041The present invention provides a 3D forward modeling method to compute gravity and tensor gravity data using Poisson's equation for airborne, surface and borehole applications in exploration problems. The forward modeling is based on a differential equation approach for solving the Poisson's (or Laplace's) equation using high order compact finite difference methods or finite element methods.

0042The formulation of the present invention makes it possible to build inversion processes to predict or determine subsurface parameters, structures and areas of interest so as to minimize the difference between any or all of the model fields and their measured counterparts. Since this inversion process is so robust and applicable to so many varied field measurements, it is useful for exploration at both regional and prospect scales.

0043U.S. Pat. Nos. 6,278,948, 6,424,918, 6,430,507, and 6,502,037 and U.S. patent application Ser. No. 10/236,204 having the same assignee and the contents of which are fully incorporated here by reference, disclose methods and apparatus in which gravity and magnetics inversion techniques developed for use with gravity, Full Tensor Gradiometry (FTG) gravity and magnetics can be used as a driver for building earth models that may be used with seismic data processing. Some of the disclosures of the '948, '918, '507 and '037 patents and the '204 application are included here for convenience.

0044These previous inventions utilize very robust formulations for calculating vector and tensor gravity and magnetic fields due to a parameterized model such as salt formations embedded in a known or assumed distribution of sediments. The embedded model is composed of an upper boundary, usually mapped from seismic data, and a parameterized lower boundary to be determined by an inversion process. The parameterized lower boundary first uses gravity and/or magnetic data to predict parameters. Then, seismic data can be combined to provide a depth image using the predicted parameters as the driver. Further, a velocity model can be derived using the predicted parameters. In the alternative, the seismic data may be used as an additional constraint in running the inversion process again. This process of inversion followed by seismic imaging followed by another inversion and seismic imaging step may be repeated until the results of the potential fields inversion and the seismic imaging processes converge to a single answer or a point of diminishing returns is reached. Otherwise, if the process results begin to diverge indicating that there is not a unique solution, the process is discontinued.

0045The inversion process of the present invention demonstrates its strength in several areas. In the area of background sediment and subsurface formation properties, the density can have any horizontal variability including discontinuities. The depth variation of these properties is assumed to be a polynomial up to order six for density. Thus, for practical purposes, there are no apparent restrictions on the variability of density especially considering the inherent resolution of gravity analyses.

0046Computations based on parallel computing techniques are very fast for three-dimensional models. The inversion process setup is ideal for using multiple processors. Thus, the elapsed time previously required for large three dimensional inversion processes can be greatly reduced. This is especially important for regional studies where large data sets are required.

0047For data preparation, the forward modeling portion of this inversion method is also used to correct the gravity and gravity tensor data for complex bathymetric effects. Bouguer corrections are made assuming a depth dependant density in order to obtain better input gravity and tensor data for the inversion process. This represents a substantial improvement over the standard use of a constant Bouguer density.

0048Some of the major problems confronting the application of any modeling technique applied to potential field data are the removal of the regional field and, thus, the isolation of the field believed to be associated with the scale of the model being examined. By properly filtering the data between successive iterations in the inversion process, predicted regional fields can be obtained for gravity and magnetic data sets simultaneously, thus allowing convergence to a common model structure that explains both the band-limited gravity and magnetic fields associated with the anomalous body. The resulting regional fields then can be modeled to predict deep-seated structures beneath the area of immediate interest.

0049Part of the inversion process that dramatically improves the convergence and efficiency of the algorithm is a predictive filtering procedure that reconstructs the regional field from the inversion itself. At each inversion step, the inversion estimates the anomalous body model that is required to fit the data and compares this model response to the observed field. The difference between the model and the observed field is treated as a “residual” or “error”, the long wavelength component of this error field is calculated and attributed to the regional field that must be accounted for in the regional model. This long wavelength residual is used to reconstruct the regional model and this reconstructed regional model is compared to the long-wavelength regional component that is removed early in the preprocessing of the potential fields data to make sure that the signal used in the inversion was properly separated between signal related to the anomalous bodies and that related to the regional field.

0050An important aspect of this invention is that the observation positions can be distributed arbitrarily, and even be positioned inside the source area, thus removing the need to have the data in a gridded form on the same datum. This is important when simultaneously inverting different data sets and using borehole gravimeter data. Also, data sets from various contractors may have very different acquisition characteristics including geographic grids. Additionally, individual field observations may be weighted in the inversion process according to their uncertainty.

0051Gravity data are generally acquired for the exploration of oil and gas or other natural resources like minerals and ore bodies. The different data types that are commonly acquired depend on the exploration objectives. The different types are: surface gravity data, tensor gravity data, airborne gravity data and borehole gravity data. The goal is to obtain subsurface density information or structural information from these different data sets by forward modeling and inversion. In principle, the density anomalies in the earth produce signals recorded in these data sets. Therefore one of the objectives is to compute the earth gravity response to subsurface anomalies using forward modeling with a known or estimated density model. The current invention of solving the Poisson's equation (using a high-order compact finite difference method) allows us to generate the gravity and tensor gravity responses in the entire medium as well as in the air with a single forward modeling run. This allows us to compute surface and borehole gravity responses concurrently without any additional complication in forward modeling. Prior art practice has been to compute the gravity and tensor gravity responses using an integral formulation that is much more complex, especially in the regions containing the source. This also requires more memory and computation time. The forward modeling approach method of the present invention provides the flexibility to invert surface, airborne and borehole gravity data simultaneously to recover density distribution in the subsurface.

0052As is well known by practitioners in the art, the divergence theorem may be applied to the gravitational field ‘g’: <maths id="MATH-US-00001" num="00001"><math overflow="scroll"><mrow><mrow><msub><mo>∫</mo><mi>V</mi></msub><mo>⁢</mo><mrow><mrow><mi>Δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>·</mo><mi>g</mi></mrow><mo>⁢</mo><mrow><mo>ⅆ</mo><mi>v</mi></mrow></mrow></mrow><mo>=</mo><mrow><msub><mo>∫</mo><mi>S</mi></msub><mo>⁢</mo><mrow><msub><mi>g</mi><mi>n</mi></msub><mo>⁢</mo><mrow><mrow><mo>ⅆ</mo><mi>S</mi></mrow><mo>.</mo></mrow></mrow></mrow></mrow></math></maths><img file="US6993433B2_D0001.tif" /><br /> With no attracting matter within the volume, ∇·g=0. This can be written as <br />∇·∇<i>U=∇</i>²<i>U=<b>0</b>.</i>

0053Prior art methods for modeling of tensor gravity data are based on solving the integral formulation for the different components of tensor gravity and a component for the vertical attraction. These methods can be used easily when the point of observation is outside the source to calculate gravity components on the surface. The methods involve the integration of each differential volume into which a model is discretized.

0054If there is a particle of mass ‘m’ within the volume ‘v’ then the surface integral is different than zero and equal to: <maths id="MATH-US-00002" num="00002"><math overflow="scroll"><mrow><mrow><msub><mo>∫</mo><mi>S</mi></msub><mo>⁢</mo><mrow><msub><mi>g</mi><mi>n</mi></msub><mo>⁢</mo><mrow><mo>ⅆ</mo><mi>S</mi></mrow></mrow></mrow><mo>=</mo><mrow><mrow><mo>-</mo><mn>4</mn></mrow><mo>⁢</mo><mi>π</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>γρ</mi></mrow></mrow></math></maths><img file="US6993433B2_D0002.tif" /><br /> This gives Poisson's equation <br />∇²<i>U=−</i>4πγρ,<br /> where ρ is the source density and γ is a constant (the gravitational constant in the case of mass and gravitational potential). Poisson's equation may be written alternatively as, <br />∇²<i>U</i>(<i>x, y, z</i>)=−4πγρ(<i>x, y, z</i>).

0055The forward modeling of Poisson's equation using the method of the present invention provides the flexibility to generate data in the entire medium including the source regions and thus allows inversion of surface and borehole gravity to recover density distribution in the subsurface. This is accomplished using the High-order compact finite difference (HOCFD) method. From a modeling standpoint this allows for generation of gravity and tensor gravity fields for the surface as well as for borehole data.

0056Rewriting Poisson's equation as <br />∇²<i>φ=ƒ,</i><br /> the equation may be discretized in this form as <br />δ_x²φ_ijk+δ_y²φ_ijk+δ_z²φ_ijk−τ_ijk<i>=ƒ</i>_ijk.<br /> The term τ_ijkis the error from the approximation. Central differences are used to <br /> calculate the operator in the discrete form using <maths id="MATH-US-00003" num="00003"><math overflow="scroll"><mrow><mrow><msubsup><mi>δ</mi><mi>x</mi><mn>2</mn></msubsup><mo>⁢</mo><msub><mi>ϕ</mi><mi>ijk</mi></msub></mrow><mo>=</mo><mrow><mfrac><mrow><msub><mi>ϕ</mi><mrow><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow><mo>,</mo><mi>j</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>-</mo><mrow><mn>2</mn><mo>⁢</mo><msub><mi>ϕ</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>,</mo><mi>k</mi></mrow></msub></mrow><mo>+</mo><msub><mi>ϕ</mi><mrow><mrow><mi>i</mi><mo>-</mo><mn>1</mn></mrow><mo>,</mo><mi>j</mi><mo>,</mo><mi>k</mi></mrow></msub></mrow><mrow><mn>2</mn><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>Δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow></mfrac><mo>.</mo></mrow></mrow></math></maths><img file="US6993433B2_D0003.tif" /><br /> Performing an expansion then: <maths id="MATH-US-00004" num="00004"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>x</mi><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mi /><mo>⁢</mo><mrow><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow><mo>+</mo><mrow><mfrac><mrow><mo>∂</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>2</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>3</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>3</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>3</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>3</mn></msup></mrow><mo>+</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi /><mo>⁢</mo><mrow><mrow><mfrac><mrow><msup><mo>∂</mo><mn>4</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>4</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>4</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>4</mn></msup></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>5</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>5</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>5</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>5</mn></msup></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>6</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>6</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>6</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>6</mn></msup></mrow></mrow></mrow></mtd></mtr></mtable></math></maths><maths id="MATH-US-00004-2" num="00004.2"><math overflow="scroll"><mi>and</mi></math></maths><maths id="MATH-US-00004-3" num="00004.3"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>x</mi><mo>-</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mi /><mo>⁢</mo><mrow><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow><mo>-</mo><mrow><mfrac><mrow><mo>∂</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>2</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup></mrow><mo>-</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>3</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>3</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>3</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>3</mn></msup></mrow><mo>+</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi /><mo>⁢</mo><mrow><mrow><mfrac><mrow><msup><mo>∂</mo><mn>4</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>4</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>4</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>4</mn></msup></mrow><mo>-</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>5</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>5</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>5</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>5</mn></msup></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>6</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>6</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>6</mn></msup></mrow></mrow></mfrac><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>6</mn></msup></mrow></mrow></mrow></mtd></mtr></mtable></math></maths><br /> Adding these equations and dividing by the cell thickness leaves <maths id="MATH-US-00005" num="00005"><math overflow="scroll"><mrow><msub><mi>ϕ</mi><mi>xx</mi></msub><mo>=</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>2</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mrow></mfrac><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>4</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>4</mn></msup></mrow></mfrac><mo>⁢</mo><mfrac><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup></mrow><mn>12</mn></mfrac></mrow><mo>+</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>6</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>6</mn></msup></mrow></mfrac><mo>⁢</mo><mfrac><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>4</mn></msup></mrow><mn>360</mn></mfrac></mrow><mo>+</mo><mrow><mi>O</mi><mo>⁡</mo><mrow><mo>(</mo><msup><mrow><mo></mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow><mo></mo></mrow><mn>6</mn></msup><mo>)</mo></mrow></mrow></mrow></mrow></math></maths><img file="US6993433B2_D0004.tif" />

0057If all the terms in the Taylor series expansion are used, a high order approximation results. Common finite-difference methods eliminate the higher order terms leaving a second order approximation as <maths id="MATH-US-00006" num="00006"><math overflow="scroll"><mrow><msub><mi>ϕ</mi><mi>xx</mi></msub><mo>=</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mrow><mn>2</mn><mo>!</mo></mrow><mo>⁢</mo><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mrow></mfrac><mo>.</mo></mrow></mrow></math></maths><img file="US6993433B2_D0005.tif" /><br /> The error approximation term then is given by <maths id="MATH-US-00007" num="00007"><math overflow="scroll"><mrow><msub><mi>τ</mi><mi>ijk</mi></msub><mo>=</mo><mrow><mfrac><msup><mi>h</mi><mn>2</mn></msup><mn>12</mn></mfrac><mo>+</mo><mrow><mo>[</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>4</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>4</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>4</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>y</mi><mn>4</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>4</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>4</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><mrow><mi>O</mi><mo>⁡</mo><mrow><mo>(</mo><msup><mrow><mo></mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>x</mi></mrow><mo></mo></mrow><mn>6</mn></msup><mo>)</mo></mrow></mrow></mrow></mrow></math></maths><img file="US6993433B2_D0006.tif" /><br /> where the value of the fourth order derivatives is given by <maths id="MATH-US-00008" num="00008"><math overflow="scroll"><mrow><mrow><mfrac><msup><mo>∂</mo><mn>2</mn></msup><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>⁢</mo><mrow><mo>(</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mfrac><msup><mo>∂</mo><mn>2</mn></msup><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>⁢</mo><mrow><mrow><mo>(</mo><mrow><mi>f</mi><mo>-</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>y</mi><mn>2</mn></msup></mrow></mfrac><mo>-</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mo>⁢</mo><mi>ϕ</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mo>.</mo></mrow></mrow></mrow></math></maths><img file="US6993433B2_D0007.tif" />

0058Substituting the error approximation into the original discretization results in the “High Order Compact Finite Difference” equation: <maths id="MATH-US-00009" num="00009"><math overflow="scroll"><mrow><mrow><mrow><mo>[</mo><mrow><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mfrac><msup><mi>h</mi><mn>2</mn></msup><mn>6</mn></mfrac><mo>⁢</mo><mrow><mo>(</mo><mrow><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mrow><mo>)</mo></mrow></mrow></mrow><mo>]</mo></mrow><mo>⁢</mo><msub><mi>ϕ</mi><mi>ijk</mi></msub></mrow><mo>=</mo><mrow><msub><mi>f</mi><mi>ijk</mi></msub><mo>+</mo><mrow><mrow><mfrac><msup><mi>h</mi><mn>2</mn></msup><mn>12</mn></mfrac><mo>⁡</mo><mrow><mo>[</mo><mrow><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mrow><mo>]</mo></mrow></mrow><mo>⁢</mo><msub><mi>f</mi><mi>ijk</mi></msub></mrow><mo>+</mo><mrow><mi>O</mi><mo>⁡</mo><mrow><mo>(</mo><msup><mi>h</mi><mn>4</mn></msup><mo>)</mo></mrow></mrow></mrow></mrow></math></maths><img file="US6993433B2_D0008.tif" /><br /> This equation corresponds to a 19-point stencil as illustrated in FIG. <b>2</b>A and FIG. <b>2</b>B. In this equation, the right hand side is made up of the usual approximation δx²+δy²+δz²and the higher order contribution <maths id="MATH-US-00010" num="00010"><math overflow="scroll"><mrow><mfrac><msup><mi>h</mi><mn>2</mn></msup><mn>6</mn></mfrac><mo>⁢</mo><mrow><mrow><mo>(</mo><mrow><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup><mo>⁢</mo><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mrow><mo>)</mo></mrow><mo>.</mo></mrow></mrow></math></maths><img file="US6993433B2_D0009.tif" /><br /> The left hand side is made up of the source function ƒ_ijkand the second derivative of the source function, <maths id="MATH-US-00011" num="00011"><math overflow="scroll"><mrow><mrow><mfrac><msup><mi>h</mi><mn>2</mn></msup><mn>12</mn></mfrac><mo>[</mo><mrow><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>x</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>+</mo><mrow><mi>δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mrow><mo>]</mo></mrow><mo>⁢</mo><mrow><msub><mi>f</mi><mi>ijk</mi></msub><mo>.</mo></mrow></mrow></math></maths><img file="US6993433B2_D0010.tif" />

0059The 19-point High Order Finite Difference Stencil “L” is illustrated in FIG. <b>2</b>A and <figref idref="DRAWINGS">FIG. 2B</figref> encompassing all the nodes of the mesh located on the three grid planes which intersect the center node of the cube, but not the corner points of the surrounding cube. An expanded view of the stencil with its matrix coefficients is illustrated in FIG. <b>2</b>A. The 5 solid vertices of the 9 point planes <b>201</b> and <b>205</b> of the stencil are combined with the 9 vertices on plane <b>203</b> to make up a 19 point finite difference stencil illustrated in FIG. <b>2</b>B. This is a 19 point operator for equi-spaced grid points. Stencils for other grids are known in the art and may be incorporated into the method of the present invention. For example, a 27 point operator formulation (for an O(h⁶) method) is described in ‘A High-Order Compact Formulation for the 3D Poisson Equation’ by W. F. Spotz and G. F. Carey, Numerical Methods for Partial Differential Equations, 12, 1996, pp. 235-243.

0060Stencils for arbitrarily spaced grids are known in the art and may be incorporated into the method of the present invention. These arbitrarily spaced grids will allow for varying boundary conditions.

0061Poisson's equation ∇²U=−4πγρ may be expressed in matrix form as <br /><i>LU</i>_ijk=−4πγρ_ijk<br /> where L is the Finite Difference diagonal dominant matrix. The finite difference matrix can be stored in a compact sparse row format that requires only three vectors, therefore reducing memory significantly over other methods. The Finite Difference system of equations can be inverted using Conjugate Gradient methods and their variants that require one matrix vector multiplication. A conjugate gradient method implementation works quickly for compact sparse row matrices compared to other methods.

0062The conjugate gradient method is a method for finding the nearest local minimum (the smallest value of a set, function, etc., within some local neighborhood) of a function of n variables which presupposes that the gradient of the function can be computed. A gradient in three dimension is, for example, the vector sum of the partial derivatives with respect to the three-coordinate variables x, y, and z of a scalar quantity whose value varies from point to point. The conjugate gradient method uses conjugate directions instead of the local gradient for going downhill. If the vicinity of the minimum has the shape of a long, narrow valley, the minimum is reached in far fewer steps than would be the case using other methods, for example the method of steepest descent.

0063As an example of an embodiment of the present invention, a discretized operator L in the expression L{right arrow over (U)}={right arrow over (f)} may be computed and examined for comparison using an analytic function for which the exact solution is known. Here <ul id="ul0003" list-style="none"><li id="ul0003-0001" num="0000"><ul id="ul0004" list-style="none"><li id="ul0004-0001" num="0064">U=sin πx sin πy sin πz,</li><li id="ul0004-0002" num="0065">ƒ=−3π²sin πx sin πsin πz, for</li><li id="ul0004-0003" num="0066">Ωε[0,1]³.</li></ul></li></ul>

0067<figref idref="DRAWINGS">FIG. 3</figref> illustrates a comparison of the analytic solution for the observed data from the potential U on the left-hand side <b>301</b>, <b>303</b> with the calculated data on the right <b>305</b>, <b>307</b>. The lower panels <b>303</b> and <b>307</b> illustrate the potential U at the center of the volume.

0068<figref idref="DRAWINGS">FIG. 4A</figref> illustrates an example of forward modeling the gravity fields using <br />∇²(Δ<i>U</i>)=−4πγΔρ

0069where the density change Δρ=−0.4gr/cm³. The integral formulation results of the model are illustrated in <figref idref="DRAWINGS">FIG. 4A</figref> with the vertical gravity g_z, and gravity tensors g_xx, g_yy, g_zz. The integral formulation results may be favorably compared with the results of the differential formulation as illustrated in <figref idref="DRAWINGS">FIG. 4B</figref> with the vertical gravity g_z, and gravity tensors g_xx,g_yy,g_zz. <figref idref="DRAWINGS">FIG. 4C</figref> illustrates the gravity tensors g_xz, g_yz, g_xyfor the integral formulation and <figref idref="DRAWINGS">FIG. 4D</figref> illustrates the gravity tensors g_xz, g_yz, g_xyfor the differential formulation of the present invention. It is important also to note the boundary condition that the potential (U) goes to zero as distance goes to infinity. <figref idref="DRAWINGS">FIG. 4E</figref> illustrates the convergence of the system of equation using conjugate gradient method. It took 10 min to solve on a SUN blade machine with ¾ million nodes.

0070<figref idref="DRAWINGS">FIG. 5A</figref> illustrates a realistic model with parameters which approximate a salt diapir with a density contrast of 0.4 grams per cubic centimeter, as shown by the change in density scale on the right. <figref idref="DRAWINGS">FIG. 5B</figref> illustrates the results of the method of the differential formulation of the present invention used to reconstruct the diapir with the vertical gravity g_z, and gravity tensors g_xx, g_yy, g_zzand g_xz, g_yz, g_xy.

0071As the method of the present invention can be used inside the volume containing the sources and inside the sources, this method has application to borehole gravity data. Prior art borehole gravity technology calculates an apparent density from the observed vertical gravity as follows: <maths id="MATH-US-00012" num="00012"><math overflow="scroll"><mrow><msub><mi>G</mi><mi>zz</mi></msub><mo>=</mo><mrow><mfrac><mrow><mrow><msub><mi>g</mi><mi>z</mi></msub><mo>⁡</mo><mrow><mo>(</mo><msub><mi>z</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mo>-</mo><mrow><msub><mi>g</mi><mi>z</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mrow><msub><mi>z</mi><mn>0</mn></msub><mo>+</mo><mrow><mi>Δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>z</mi></mrow></mrow><mo>)</mo></mrow></mrow></mrow><mrow><mi>Δ</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>z</mi></mrow></mfrac><mo>.</mo></mrow></mrow></math></maths><img file="US6993433B2_D0011.tif" /><br /><figref idref="DRAWINGS">FIG. 6</figref> illustrates this Δz dependence schematically. Apparent density is calculated as: <maths id="MATH-US-00013" num="00013"><math overflow="scroll"><mrow><msubsup><mi>ρ</mi><mi>a</mi><mi>pred</mi></msubsup><mo>=</mo><mrow><mfrac><mrow><mo>-</mo><mn>1</mn></mrow><mrow><mi>r</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>π</mi><mo>⁢</mo><mstyle><mtext> </mtext></mstyle><mo>⁢</mo><mi>γ</mi></mrow></mfrac><mo>⁢</mo><mrow><msub><mi>G</mi><mi>zz</mi></msub><mo>.</mo></mrow></mrow></mrow></math></maths><img file="US6993433B2_D0012.tif" />

0072As illustrated in <figref idref="DRAWINGS">FIG. 7</figref>, calculation of apparent density from borehole gravity using only the G_zzcomponent <b>702</b> of the field results in an anomalously high density value AV. Combining the effect due to G_xxand G_yyas illustrated in panel <b>703</b> results in a more accurate value for the density anomaly as shown in panel <b>701</b>. The contribution of each tensor component to the apparent density is illustrated in FIG. <b>8</b>. Panel <b>801</b> illustrates component G_xxnormalized by G_xx+G_yy+G_zz. Panel <b>803</b> illustrates component G_yynormalized by G_xx+G_yy+G_zz. Panel <b>805</b> illustrates component G_zznormalized by G_xx+G_yy+G_zz.

0073When the tensor fields are modeled with the differential form of the present invention we get, as a bonus, the capability of doing upward continuations as part of the calculations. The upward continued field is a result from computing the forward model.

0074In a preferred embodiment of the present invention the forward modeling method for tensor gravity data is implemented by solving Poisson's equation using a 4^thorder accuracy finite difference scheme. This enables accurate computation of the tensor gravity components. This implementation has the power of modeling surface tensor gravity data as well as borehole gravity data. Using horizontal tensor components of borehole gravity data contributes to a better prediction of the apparent density. Accurate modeling of borehole gravity data allows for more efficient monitoring of oilfields over time. The method adapts easily to compute upward and downward continued potential fields.

0075Other variations are within the scope of the present disclosure. For example, the method is readily adaptable to be used within an inversion procedure that uses integral equations to calculate a forward model.

0076Gravity, tensor gravity and borehole data are acquired in the exploration for oil, gas and minerals. These data provide a means to image the subsurface, based on density anomalies in the earth. The high order compact finite difference methodology using this work allows us to obtain accurate gravity (third order accuracy i.e. O(h³)) and tensor gravity (second order accuracy, i.e. O(h²)) in the entire medium. This methodology significantly improves the accuracy over the conventional nine point finite difference approaches and therefore is highly advantageous for the modeling these field. The finite difference stencil used in a preferred embodiment is a 19 point stencil. The resulting system of equations is stored in a compressed row sparse format for use with the conjugate gradient inversion techniques as described earlier.

0077The gravity (G_x, G_y, G_z) and gravity tensor components (G_xx, G_yy, G_zz, G_xy, G_xz, G_yz) are computed by obtaining the gravitational potential in the entire discretized medium using the above method. Subsequently the gravity and the tensor components are evaluated using forward and central difference formulas. The values of the components that are staggered with respect to the original mesh (coordinate layout) are interpolated to the initial nodal coordinates using a 3D interpolation method. The method of the present invention also provides for a formulation that can handle unequal cell sizes to compute the potential. Adapting the present High Order Compact formulation for use with unequal cell sizes has the advantage of imposing the Dirichlet boundary condition (i.e. where the potential goes to infinity) with fewer cells.

0078Referring now to <figref idref="DRAWINGS">FIG. 9</figref>, this flow chart illustrates a method for providing an initial model upon which derivation of subsurface formations (the geophysical model) can be based. At block <b>12</b> reflection seismic data is received. In general, reflection seismic data is more reliable than other forms of geologic exploration data in most situations because of the extensive work and research that has been devoted to its use and application.

0079At block <b>14</b>, the reflection seismic data is used to derive the top portion of a geologic model. As stated previously, reflection seismic data is very reliable, at least in this situation, for determining a model containing the top of an anomalous density zones. Geologic boundaries beneath the top of an anomalous zone or boundary are not as easily modeled. While reflection seismic data is generally more reliable than other forms of geological surveying, anomalous density zones, such as a salt dome, obstruct, divert or highly distort seismic reflection information below the boundary or zone. In this type of area, reflection seismic surveying becomes unreliable, as shown in the example of FIG. <b>1</b>. In some cases, the reflection seismic data becomes so unreliable that no useful information below the salt dome can be obtained from reflection seismic exploration.

0080At block <b>16</b> data pertaining to the determination of the lower boundary is received. This data may take the form any potential fields data, both vector and tensor. In the formulation used herein, the type of data is of no concern since any combination of the above mentioned data can be processed simultaneously. Although these types of data generally provide less resolution than reflection seismic data, in the case of an anomalous density zone, potential field data may provide the most reliable available data.

0081The data received are used to formulate the limits of the lower boundary for the geologic model block <b>18</b>. The actual lower boundary derivation is done by predicting parameters representing the lower boundary, and this predicting maybe done by a-priori knowledge or through an inversion process.

0082Although various inversion techniques can be utilized to determine the parameters (coefficients) representing the lower boundary, one preferred embodiment of the present invention involves the successive inversion of a single coefficient at a time until all coefficients are determined. The total number of coefficients is set apriori at block <b>18</b> based upon the minimum wavelength (maximum frequency) desired in the lower boundary. Typical lower boundaries may contain as many as nine hundred coefficients for three-dimensional models, thirty in each horizontal direction. For example, if x₁and x₂are the spatial limits of integration in the x direction, and y₁and y₂are the spatial limits of integration in the y direction, and half cosine series are used, the number of terms required in the x and y directions respectively for a minimum wavelength of λ_minare: <maths id="MATH-US-00014" num="00014"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mi>n</mi><mi>x</mi></msub><mo>=</mo><mrow><mfrac><mrow><mn>2</mn><mo>⁢</mo><mrow><mo>(</mo><mrow><msub><mi>x</mi><mn>2</mn></msub><mo>-</mo><msub><mi>x</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow></mrow><msub><mi>λ</mi><mi>min</mi></msub></mfrac><mo>+</mo><mn>1</mn></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mi>and</mi></mtd><mtd><mstyle><mtext> </mtext></mstyle></mtd></mtr><mtr><mtd><mrow><msub><mi>n</mi><mi>y</mi></msub><mo>=</mo><mrow><mfrac><mrow><mn>2</mn><mo>⁢</mo><mrow><mo>(</mo><mrow><msub><mi>y</mi><mn>2</mn></msub><mo>-</mo><msub><mi>y</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow></mrow><msub><mi>λ</mi><mi>min</mi></msub></mfrac><mo>+</mo><mn>1</mn></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><img file="US6993433B2_D0013.tif" /><br /> Thus the total number of coefficients representing the lower boundary is n_xmultiplied n_y.

0083At block <b>20</b>, the first coefficient representing a uniform lower boundary (longest wavelength component) is predicted. This coefficient is based on the received data and represents the best fit possible with limited information.

0084At block <b>22</b> another coefficient is added through the inversion process. At block <b>24</b> all coefficients added thus far are readjusted through single parameter inversion in the order in which they were added to the lower boundary (from the longest to the shortest wavelength coefficients). This process continues until all coefficients (n_xtimes n_y) have been predicted. At decision block <b>26</b> a determination is made as to whether all coefficients have been predicted. If they have not, the program returns to block <b>22</b> where another coefficient is added.

0085If all coefficients have been predicted, the program proceeds to block <b>28</b> where the lower boundary is displayed and subsequently sent on to the forward modeling of the present invention. The display may take any form currently in use in the art, such as a monitor display, a hard copy printout or be kept electronically on tape or disc.

0086Initial geophysical model development and the updating of geophysical models based the differences between measured and estimated gravity and gravity tensor data are extensively covered in the disclosures of the '948, '918, '507 and '037 patents and the '204 application which have been fully incorporated herein by reference. Further updating geophysical models based on integration of potential fields and seismic data is covered in these disclosures as well.

0087A number of embodiments of the '204 application are described therein which are applicable for combination with method of the present invention. In one embodiment, the predicted parameters are combined with seismic data to obtain a depth image and a velocity model for delineating formations or parameters of interest. In a second embodiment, imaging of the seismic data is performed using a lower boundary determined from the initial inversion. This seismic result is mapped to determine the position of the base of the anomalous body and is compared to the model parameters to obtain a difference between the two. The parameters are adjusted to provide a best fit. A new model is formed based on the new predicted parameters, and the potential fields modeling (using HOCFD in the present case) and seismic imaging modeling results are iterated until the solution converges.

0088Where an anomalous body or subterranean formation of interest is a volcanic intrusive or extrusive body that has higher density than the surrounding sedimentary formations, the overburden stress will be higher than would be expected in a comparable depth of sedimentary rocks. The present invention is able to account for the effects of zones of anomalous density by determining the actual density, whereas prior art methods assume that density increases uniformly with depth or varies with depth based on limited well information.

0089One embodiment of the present invention includes a modification of the very fast inversion process which produces models using vector and tensor potential methods (applicable for both gravity and magnetic data). This modification involves solving the inverse problem by using a constrained optimization approach, for example, a constrained nonlinear inversion method. The formulation allows imposing constraints on the solution as well as incorporating any geologic apriori information in the inversion algorithm. The computational speed is achieved using a conjugate gradient solver and using approximate sensitivity when required. The algorithm has the flexibility to allow an interpreter to carry out hypothesis testing, based on his geologic intuition or any other information about the model derived from any sources, for example, geophysical or geological data.

0090The method of this embodiment applies a constrained nonlinear inversion method using a Gauss-Newton approach. The model is parameterized using, for example, a boxcar basis function. Then the forward problem is solved using the method described above. In the inversion algorithm we invert for all model parameters simultaneously and optimally constrain parameters by imposing soft bounds on the solution. This bounded minimization problem is solved using an interior-point method. Interior-point methods are well known in linear programming. The interior-point algorithm guarantees that the solution will remain within bounds, i.e. the inverted base of salt will always be greater than (or equal to) the top of salt and less than any prescribed upper bound. This formulation has many model parameters compared to the data, which leads to an underdetermined problem. Therefore, to address the non-uniqueness issues, we minimize a model objective function subject to fitting the data. The choice of this model objective function is generally geologically driven. For example, we can impose smoothness criteria in some regions of the model and also allow sharp features to be built in other regions like faults that may be based on knowledge or geological objectives. In essence this model objective function can be used to impose a constraint on the model based on apriori knowledge. For example, a topological constraint may be placed on the base of the salt model (or any other parameter) based on apriori knowledge, and geological or geophysical information of the subsurface.

0091Since we are solving for all of the model parameters in this approach, the system of equations may require a large matrix system be solved. Using a least-squares conjugate gradient solver and using approximate sensitivity when necessary achieves an increase in computation speed. The approximate sensitivity requires less storage and allows the matrix to be stored in a sparse form with very minimal loss in the information content due to approximation. The noise in the data is handled using data weighting matrices. Appropriate stopping criteria for the convergence of the inversion process may be implemented and based on the assumption that the noise is random and has a Gaussian distribution. Other distributions and parameterization functions may be easily incorporated into this methodology.

0092The constrained nonlinear inversion method using a Gauss-Newton approach leads naturally to an addressing of resolution issues for the inverted model. In the regions that are not constrained by the data, the model parameters from inversion may not be realistic. To avoid this problem, a reference model may be built into the inversion algorithm. Thus if the regions are insensitive to the data, the recovered model from inversion will reflect the reference model. This can be used to identify regions in the inverted model that are demanded by the data and thus provide a better interpretability of the results.

0093The method of this embodiment is further explained with reference to the flow chart of FIG. <b>10</b>. Gravity data <b>401</b>, which may include gravity tensor data, are acquired for input to the method of the invention. An initial geophysical model <b>403</b> is determined or derived from apriori knowledge and/or other available data. A model potential field <b>405</b> is determined from the initial geophysical model <b>403</b>. The initial geophysical model can be made using the forward modeling of the High Order Compact Finite Difference (HOCFD) method of the present invention. A difference is determined <b>407</b> between the acquired gravity data <b>401</b> and the gravity data potential fields' response due to the geophysical model <b>405</b>. If the difference <b>407</b> is outside of a selected minimum range <b>409</b>, an inversion process <b>410</b> is implemented, for example through applying a constrained nonlinear inversion using a Gauss-Newton approach as outlined in the '204 application. The model <b>405</b> is parameterized using appropriate basis functions; model parameters are inverted for simultaneously; soft bounds are imposed on the solution; the minimization problem may be solved with an interior-point method. The non-uniqueness of the underdetermined problem is addressed by minimizing one or more model objective functions, for example, by imposing geological constraints based on apriori knowledge or interpretations. In this manner, the estimated gravity potential fields' response <b>405</b> is compared with measured gravity data <b>401</b> and iterated through <b>410</b> until the model differences <b>407</b> are within selected limits <b>409</b>.

0094Two other embodiments of the present invention have application to address prior art data noise problems with the use of gravity tensor technology. One of the major problems with tensor methods is the typically high noise level associated with most tensor surveys. In fact, in many cases the noise may swamp the signals believed to be attributable to geologic structures. Before the tensor data can be used for exploration purposes, a process must be applied to try to extract the signal from the measured tensor data.

0095Embodiments of the present invention for application to noise problems may incorporate one or both of two techniques that may be applied separately or in tandem to drastically reduce the level of noise in all five tensor channels. These two techniques are the ‘Laplace's Equation’ method and the ‘equivalent source’ method. Barely interpretable datasets can be converted to datasets that may be almost ‘geologic’ in nature.

0096The method based on a solution of Laplace's equation works very well for merging the low spatial frequencies measured with a standard gravimeter with the entire spatial frequency spectrum measured with full gravity tensors measurements. The Laplace's equation method utilizes a solution on an arbitrary surface not intersecting the region containing the sources causing the gravity and tensor response. As an example, the method may use a 2D half-cosine transform-like function for the x and y dependence, along with the exponential function for the z dependence. Thus, the gravitation potential V(x,y,z), can be represented as a series of products of the two functions with coefficients A_mnwhere m and n refer to the spatial wave numbers Kx, Ky and Kz in the x, y and z direction respectively. Kz becomes a function of Kx and Ky when this series is made to be a solution to Laplace's equation for V(x,y,z). The utility of this form of V(x,y,z) is that it can be easily included in a linear-inverse scheme to investigate the true signal in a set of gravity and gravity tensor measurements. Since the gravity field and tensor can be derived from the gravitational potential V(x,y,z), the A_mn's can thus be estimated by fitting the observed data using various linear-inversion algorithms. This allows the components of all data channels that satisfy Laplace's equation simultaneously to be isolated from the noise. Noise may be defined as that part of the measurements that does not pass this test. The entire survey can be condensed down to a single set of coefficients A_mn. These coefficients can be used for producing arbitrary grids of all channels. In addition, various operations can be performed such as upward and downward continuation, higher-order derivative maps, etc. Furthermore, 3D modeling such as that outlined previously can be greatly simplified using output from this process. A case in point is that the modeling program can use one channel such as the second vertical derivative tensor (G_zz), since it is directly calculated from V(x,y,z), and therefore contains information from all field measurements. This fact drastically reduces CPU time and memory requirements for doing 3D modeling. Hence, the method provides the opportunity of applying more sophisticated modeling techniques to larger survey areas while using fewer resources. One further point is that this method can be applied to randomly distributed observation locations since no FFT operations are used and the whole process is model independent.

0097The Laplace's Equation method embodiment is further explained with reference to the flow chart of FIG. <b>11</b>. Potential fields data <b>401</b> are acquired. An initial geophysical model <b>403</b> is determined or derived from apriori knowledge and/or other available data. A model potential field as a solution to Laplace's Equation method <b>405</b> is determined from the initial geophysical model <b>403</b>, using for example, half-cosine function transforms in x and y, and an z variable exponential function. A difference is determined <b>407</b> between the acquired potential fields data <b>401</b> and the potential fields response due to the geophysical model <b>405</b>. If the difference <b>407</b> is outside of a selected minimum range <b>409</b>, a forward modeling process <b>420</b> is implemented by applying a Laplace's Equation solution approach and/or an Equivalent Source method. In this manner, the estimated or modeled potential fields response <b>405</b> is compared with measured potential fields data <b>401</b> and iterated <b>410</b> until the model differences <b>407</b> are within selected limits <b>409</b>.

0098Another method that may be used for data noise reduction is based on equivalent-source modeling. Equivalent-source modeling deals with the formation of an infinitely thin layer of spatially varying mass between the source(s) believed to be responsible for the measured gravity and gravity tensor fields, and the observation positions. In the case of marine surveys this layer of surface mass may be chosen to be coincident with the bathymetric surface. The method of forward modeling has been outlined previously. The main difference is that the surface density is set as the product of a two-dimensional surface density distribution function, Ms(x′,y′), and a delta function with the argument (z′-Z(x′,y′)), where Z(x′,y′) is the bathymetry. This results in an analytical integration of z′, thus leaving a numerical 2D integration over x′ and y′ for the model field calculations. In a linear inversion program, Ms(x′,y′) can be determined that best models the measured tensors and gravity simultaneously. Ms(x′,y′) is parameterized analogously to the base of salt as previously disclosed. Just as for the Laplace equation solution method, the residual fields that cannot be modeled are considered noise, and this method also can utilize randomly distributed observation points. The Laplace equation solution method and the equivalent source modeling method may be used in tandem.

0099The Equivalent Source method is further explained with reference to the flow chart of FIG. <b>11</b>. Potential fields data <b>401</b> are acquired. An initial geophysical model <b>403</b> is determined or derived from apriori knowledge and/or other available data. A model potential field as a solution to the Equivalent Source method <b>405</b> is determined from the initial geophysical model <b>403</b>, using for example, surface density set as the product of a two-dimensional surface density distribution function, Ms(x′,y′), and a delta function with the argument (z′-Z(x′,y′)), where Z(x′,y′) is the bathymetry. A difference is then determined <b>407</b> between the acquired potential fields data <b>401</b> and the potential fields response determined from the geophysical model <b>405</b>. If the difference <b>407</b> is outside of a selected minimum range <b>409</b>, a forward modeling process <b>420</b> is implemented by applying an Equivalent Source method. In this manner, the estimated or modeled potential fields response <b>405</b> is compared with measured potential fields data <b>401</b> and iterated <b>410</b> until the model differences <b>407</b> are within selected limits <b>409</b>.

0100A preferred embodiment of the method of the present invention is further illustrated with reference to the flow chart of FIG. <b>12</b>. Potential fields data <b>401</b> are acquired. An initial geophysical model <b>403</b> is determined or derived from apriori knowledge and/or other available data using HOCFD. A model gravity data potential field as a solution to the Poisson equation 405 is determined from the initial geophysical model <b>403</b>, using for example, data values derived a locations equivalent to the potential fields' data <b>401</b>. A difference is then determined <b>407</b> between the acquired potential fields data <b>401</b> and the potential fields response determined from the geophysical model using HOCFD <b>405</b>. If the difference <b>407</b> is outside of a selected minimum range <b>409</b>, a forward modeling process <b>430</b> is implemented by applying any inversion method wherein the modeled gravity and/or gravity tensor data are inverted to obtain values consistent with the structures of the initial geophysical model <b>401</b>, or the structures and their density distributions may be altered as outlined in embodiments above. In this manner, the estimated or modeled potential fields response <b>405</b> is compared with measured potential fields data <b>401</b> and iterated <b>410</b> until the model differences <b>407</b> are within selected limits <b>409</b>.

0101While there has been illustrated and described particular embodiments of the present invention, it will be appreciated that numerous changes and modifications will occur to those skilled in the art, and it is intended in the appended claims to cover all those changes and modifications which fall within the true spirit and scope of the present invention.