Method of terrain correction for potential field geophysical survey data

Abstract

A method for terrain correction of potential field geophysical survey data measured above an examined medium having density and/or magnetization is described, using potential field data including but not limited to gravity and/or magnetic total field and/or vector and/or tensor data. The potential field sensors may measure the gravity and/or magnetic total field and/or vector and/or tensor data at least one receiving position with respect to the examined medium. The terrain of the examined medium may be described by a spatially variable analytic function of the material properties of the examined medium. The terrain response for at least one component of the measured potential field in at least one receiver location (potential field data) may be calculated using special form of surface integral over the terrain based on 3D analog of the Cauchy-type integrals. This surface integration ensures accurate representation of the terrain response.

Claims

1. A method for terrain correction, the method comprising: a. measuring at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data with at least one potential field sensor in at least one measurement position along at least one survey line; b. measuring a terrain as a digital elevation model; c. selecting a physical property function to describe a physical property distribution within the measured terrain; d. using at least one processor, forward modeling a terrain response from the digital elevation model with the physical property function for the measured at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data using a 3D Cauchy-type surface integral over a surface of the terrain; e. using the at least one processor, filtering the forward modeled terrain response to match a system response and a data processing response of the at least one potential field sensor; f. using the at least one processor, calculating the filtered forward modeled terrain response for the measured at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data; g. using the at least one processor, calculating a terrain corrected potential field data using the measured at least one component and the calculated filtered forward modeled terrain response; and h. using the terrain corrected potential field data for geophysical exploration of one or more geological areas.

2. The method of claim 1, wherein a physical property of the physical property function comprises one of density, susceptibility, and/or magnetization, representing the physical properties of the terrain containing the geological areas.

3. The method of claim 1, wherein the physical property function describing the physical property distribution within the terrain is an analytic function of space.

4. The method of claim 1, wherein the physical property function describing the physical property distribution within the terrain is constant.

5. The method of claim 1, wherein the physical property function describing the physical property distribution within the terrain may vary vertically according to a predetermined law.

6. The method of claim 1, wherein the at least one potential field sensor comprises a plurality of potential field sensors arranged in an array.

7. The method of claim 1, wherein an examined medium contains the geological areas.

8. The method of claim 1, further comprising: placing at least one potential field sensor in at least one measurement position in a survey.

9. A physical non-transitory computer readable medium having stored thereon computer executable instructions that when executed by at least one processor cause a computing system to perform a method for terrain correction, the method comprising: a. measuring at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data with at least one potential field sensor in at least one measurement position along at least one survey line; b. measuring a terrain as a digital elevation model; c. selecting a physical property function to describe a physical property distribution within the measured terrain; d. forward modeling a terrain response from the digital elevation model with the physical property function for the measured at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data using 3D Cauchy-type surface integral over a surface of the terrain; e. filtering the forward modeled terrain response to match a system response and a data processing response of the at least one potential field sensor; f. calculating the filtered forward modeled terrain response for the measured at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data; g. calculating a terrain corrected potential field data using the measured at least one component and the calculated filtered forward modeled terrain response; and h. using the terrain corrected potential field data for geophysical exploration of one or more geological areas.

10. The method of claim 9, further comprising: placing at least one potential field sensor in at least one measurement position in a survey.

11. A system for terrain correction comprising: at least one potential field sensor; and a computing system, the computing system comprising: at least one processor; and one or more physical non-transitory computer readable media having computer executable instructions stored thereon that when executed by the at least one processor, cause the computing system to perform the following: measure at least one component of the following: gravity scalar, gravity vector, gravity tensor magnetic scalar, magnetic vector, and magnetic tensor data with at least one potential field sensor in at least one measurement position along at least one survey line; measure a terrain as a digital elevation model; select a physical property function to describe a physical property distribution within the measured terrain; forward model a terrain response from the digital elevation model with the physical property function for the measured at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data using 3D Cauchy-type surface integral over a surface of the terrain; filter the forward modeled terrain response to match a system response and a data processing response of the at least one potential field sensor; calculate the filtered forward modeled terrain response for the measured at least one component of the following: gravity scalar, gravity vector, gravity tensor, magnetic scalar, magnetic vector, and magnetic tensor data; calculate a terrain corrected potential field data using the measured at least one component and the calculated filtered forward modeled terrain response; and use the terrain corrected potential field data for geophysical exploration of one or more geological areas.

12. The system of claim 11, wherein the physical properties comprise one of density, susceptibility, and/or magnetization, representing the physical properties of the terrain containing geological areas.

13. The system of claim 11, wherein the physical property function describing the physical property distribution within the terrain is an analytic function of space.

14. The system of claim 11, wherein the physical property function describing the physical property distribution within the terrain is constant.

15. The system of claim 11, wherein the physical property function describing the physical property distribution within the terrain may vary vertically according to a predetermined law.

16. The system of claim 11, wherein the at least one potential field sensor comprises a plurality of potential field sensors arranged in an array.

17. The system of claim 11, wherein an examined medium contains a geological area.

18. The system of claim 11, wherein the at least one potential field sensor is placed in at least one measurement position in a survey.

19. The system of claim 11 wherein the at least one potential field sensor comprises one or more of a gravimeter, a gravity gradiometer, a global positioning system (GPS), Laser Imaging Detection and Ranging (LIDAR) altimetry, radar altimetry, magnetometers, magnetic gradiometers, natural source electromagnetic systems, controlled-source electromagnetic systems, and spectrometers.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Exemplary embodiments of the invention will become more fully apparent from the following description and appended claims, taken in conjunction with the accompanying drawings. Understanding that these drawings depict only exemplary embodiments and are, therefore, not to be considered limiting of the invention's scope, the exemplary embodiments of the invention will be described with additional specificity and detail through use of the accompanying drawings in which:

(2) FIG. 1 illustrates an embodiment of a system for terrain correction of gravity and gravity gradiometry data including gravimeters and gravity gradiometers attached to a fixed wing aircraft that is moving along a survey line L(t) above the surface of the examined medium.

(3) FIG. 2 illustrates an embodiment of a system for terrain correction of gravity and gravity gradiometry data including gravimeters and gravity gradiometers attached to a vessel moving at some elevation along a survey line L(t) above the surface of the examined medium.

(4) FIG. 3 illustrates an embodiment of a system for terrain correction of gravity and gravity gradiometry data from a fixed wing aircraft.

(5) FIG. 4 illustrates an embodiment of a processor for terrain correction of gravity and gravity gradiometry data.

(6) FIG. 5 illustrates an embodiment of a method for terrain correction of gravity and gravity gradiometry data.

(7) FIG. 6 illustrates the 3D earth model used for calculating the terrain response of the measured gravity and gravity gradiometry data.

(8) FIG. 7 illustrates a discrete cell on the surface of the terrain used for the discrete forward solution of the 3D Cauchy-type integral equations for the gravity and gravity gradient fields.

(9) FIG. 8 illustrates the plan view of the anomalous bodies in a synthetic density model used to simulate gravity gradiometry data.

(10) FIG. 9 illustrates the plan view of the LiDAR digital elevation model over the synthetic survey area.

(11) FIG. 10 illustrates the plan view of the LiDAR digital elevation model merged with the SRTM digital elevation model over the synthetic survey area used to simulate gravity gradiometry data.

(12) FIG. 11 illustrates the free-air vertical gravity gradient (Gzz) data simulated with 1 E/Hz noise added for the synthetic density model and digital elevation model.

(13) FIG. 12 illustrates the terrain corrected vertical gravity gradient (Gzz) data recovered from the synthetic free air vertical gravity gradient (Gzz) data and merged LiDAR and SRTM digital elevation models.

DETAILED DESCRIPTION

(14) Potential field geophysical surveys may encompass marine, borehole, ground and airborne potential field measurements from moving platforms such as but not limited to airplanes, helicopters, airships, unattended aerial vehicles (UAV), borehole logging instruments, vessels, submarines, and other vehicles.

(15) Potential field geophysical data usually includes at least one component of the gravity and/or magnetic vector and/or tensor. Potential field vector data may be directly measured as either at least one vector component (e.g., vertical gravity component as measured by a gravimeter) or as the magnitude of the total field (e.g., total magnetic intensity as measured by an optically pumped magnetometer), or indirectly calculated from potential field tensor data. Potential field tensor data may be directly measured (e.g., magnetics tensors as measured by a SQUID magnetic gradiometer) or indirectly calculated from potential field vector data. Potential field tensor data may also be approximated by directly and simultaneously measuring at least one component of the potential field using at least two spatially separated instruments measuring at least one vector component of the potential field (e.g., magnetic tensor as approximated by the gradient computed between total magnetic intensities as measured by two spatially separated optically pumped magnetometers).

(16) Measured potential field data contains the linear superposition of responses from both the terrain and the subsurface geological structures located beneath it.

(17) For the purpose of interpreting potential field data, the terrain response is calculated using special form of surface integral over the surface of the elevation model based on 3D analog of the Cauchy-type integrals. This surface integration ensures accurate representation of the terrain response.

(18) The computed terrain response is filtered to match the potential field instrument's system and data processing response. For example, the system and data processing response of airborne gravity gradiometers can be effectively emulated with a sixth or seventh order low-pass Butterworth filter of length between 100 m and 1000 m.

(19) The system filtered terrain response is subtracted from the measured potential field data. For example, for gravity gradiometry, the system filtered terrain response should be calculated and subtracted from the measured free-air gravity gradiometry data.

(20) In at least one embodiment of a method disclosed herein, the terrain correction can be applied to gravity vector data measured with a gravimeter, or calculated from data measured with a gravity gradiometer.

(21) In at least one embodiment of a method disclosed herein, the terrain correction can be applied to gravity tensor data measured with a gravity gradiometer, or calculated from data measured with a gravimeter.

(22) In at least one embodiment of a method disclosed herein, the terrain correction can be applied to magnetic vector data, including total magnetic intensity, measured with a magnetometer, or calculated from data measured with a magnetic gradiometer.

(23) In at least one embodiment of a method disclosed herein, the terrain correction can be applied to magnetic tensor data measured using a magnetic gradiometer, or calculated from data measured with a magnetometer.

(24) At least one embodiment of a method disclosed herein, can be applied for subsurface imaging of geological and/or man-made structures for mineral, hydrocarbon, geothermal and groundwater exploration, in-situ mining, hydrocarbon, geothermal and groundwater resource monitoring, and tunnel and UGF detection, using potential field geophysical data measured from at least one moving platform such as but not limited to airplanes, helicopters, airships, UAV, borehole logging instruments, vessels, submarines, and vehicles.

(25) At least one embodiment of this method can be used in geophysical exploration for mineral, hydrocarbon, geothermal, and groundwater resources.

(26) At least one embodiment of this method can be used in geophysical monitoring for in-situ mining, hydrocarbon, geothermal, and groundwater resources.

(27) At least one embodiment of this method can be used in geophysical monitoring for earth systems.

(28) At least one embodiment of this method can be used for tunnels and UGF detection.

(29) One embodiment of a system for terrain correction of gravity and gravity gradiometry data is illustrated in FIG. 1. A data acquisition system 1, located on a fixed wing aircraft 2, may include gravimeters and/or gravity gradiometers 3 attached to the fixed wing aircraft that is moving along a survey line L(t) 4 at some elevation above the surface of an examined medium 5 that may be characterized by a homogeneous terrain density 6 with an anomalous density distribution 7 superimposed upon it.

(30) Another embodiment of a system for terrain correction of gravity and gravity gradiometry data is illustrated in FIG. 2. A data acquisition system 8, located on a vessel 9, may include gravimeters and/or gravity gradiometers 10 attached to the vessel that is moving along a survey line L(t) 11 at some elevation above the surface of an examined medium 12 that may be characterized by a terrain density that varies as a function of depth 13 with an anomalous density distribution 14 superimposed upon it.

(31) An embodiment of a system for terrain correction of gravity and/or gravity gradiometry data from a fixed wing aircraft is illustrated in FIG. 3. The real time volume imaging system may include general aircraft instrumentation 15, gravity system (e.g., gravimeter) 16, gravity gradient system (e.g., gravity gradiometer) 17, global positioning system (GPS) 18, LIDAR altimetry 19, radar altimetry 20, other geophysical sensor systems including but not limited to magnetometers, magnetic gradiometers, natural source electromagnetic systems, controlled-source electromagnetic systems, and spectrometers 21, data acquisition system 22, communications system 23, and processor 24, which collectively can produce free air gravity and/or gravity gradiometry data 25 and terrain corrected gravity and/or gravity gradiometry data 26.

(32) An embodiment of a processor 24 is illustrated in FIG. 4. The processor 24 may include, for example, a data and code memory 27 for storing free-air gravity and/or gravity gradiometry data from the data acquisition system 22 via the communications system 23, terrain correction computer software and terrain corrected gravity and/or gravity gradiometry data, a central processing unit 28 for executing the terrain correction computer software on the free-air gravity and/or gravity gradiometry data to generate terrain corrected gravity and/or gravity gradiometry data, a graphical user interface (GUI) 29 for displaying the free-air and/or terrain corrected gravity and/or gravity gradiometry data, and a communications system 30 for system interoperability. The processor 24 may comprise of a single processing unit or can be distributed across one or more processing units in one or more locations. The communications system 23 can include I/O interfaces for exchanging information with one or more external devices. The data and code memory 27 may comprise of a single memory device or can be distributed across one or more memory devices in one or more locations connected via the communications system 23.

(33) An embodiment of a method for terrain correction of gravity and/or gravity gradiometry data is schematically shown in FIG. 5. The processor 24 may merge high resolution terrain data such as that obtained from LiDAR measurements during a survey 31 and low resolution terrain data, such as that obtained from a public geospatial database such as SRTM 32 to form a digital elevation model 33 of the survey area and an area 10 km (or other area size) from the perimeter of the survey area. The processor 24 may also read the gravimeter's and/or gravity gradiometer's geospatial data 34 that is processed from the aircraft's navigation instrumentation. A terrain density function 35 is selected by the user, and may be expressed by any analytic function. A 3D forward processor 36 may compute the predicted terrain response. The gravimeter's and/or gravity gradiometer's system responses are emulated as a low-pass filter 37 whose parameters are chosen by the user. The processor 24 may low pass filter the predicted terrain response to produce a filtered terrain response 38 that would be measured by the gravimeters and/or gravity gradiometers. The processor 24 may read the free-air gravity and/or gravity gradiometry data 39 and subtract the filtered terrain response 38 to produce terrain corrected gravity and/or gravity gradiometry data 40.

(34) The 3D forward processor 36 computes the gravity and/or gravity gradiometry responses using Cauchy-type integral representations of the gravity and/or gravity gradient fields for the density contact surface.

(35) The concept of a three-dimensional Cauchy-type integral was introduced by Zhdanov, 1988, and is represented by the following expression:

(36) $\begin{array}{cc}{C}^S\left(r,\right)=-\frac{1}{4}{}_S\left[\left(n.Math.\right)\frac{1}{.Math.r-r^{}.Math.}+\left(n\right)\frac{1}{.Math.r-r^{}.Math.}\right]ds,& \left(1\right)\end{array}$
where S is some closed surface, =(r) is some vector function specified on S and is continuous on S, and n is a unit vector of the normal to S directed outside the domain D, bounded by the surface S. Function is called a vector density of the Cauchy-type integral C.sup.S(r, ).

(37) The Cauchy-type integral C.sup.S(r, ) can be represented in matrix notation. The vectors C.sup.S, , n, and

(38) $\frac{1}{.Math.r-r^{}.Math.}$
are presented in some Cartesian basis {d.sub.x, d.sub.y, d.sub.z}, where the z axis is directed upward, as:

(39) $C^S=C_{}^S{d}_{}$$={}_{}{d}_{}$$n=n_{}{d}_{}$$\frac{1}{.Math.r-r^{}.Math.}=-\frac{r-r^{}}{{.Math.r-r^{}.Math.}^3}=-\frac{r_{}-r_{}^{}}{{.Math.r-r^{}.Math.}^3}{d}_{}$
where r.sub.= and , , , =x, y, z, and where we use an agreement on summation that the twice recurring index indicates the summation over this index. Using these notations, equation (1) can be re-written as:

(40) $\begin{array}{cc}{C}^S\left(r,\right)=-\frac{1}{4}{}_S{}_{}{}_{}\frac{r_{}-r_{}^{}}{{.Math.r-r^{}.Math.}^3}{n}_{}ds,& \left(2\right)\end{array}$
where the four-index -symbol is expressed in terms of the symmetric Kronecker symbol .sub.:

(41) ${}_{}={}_{}{}_{}+{}_{}{}_{}-{}_{}{}_{},{}_{}=\{\begin{array}{cc}1,& =,\\ 0,& .\end{array}$

(42) As illustrated in FIG. 6, we can consider a model with the density contrast at some surface F 41, representing the topography of the earth's surface, i.e., the terrain. We can consider a domain D 42 that is infinitely extended in the horizontal directions bounded by the surface F, described by equation z=h(x,y)H.sub.0, and a horizontal plane P 43, z=h.sub.0, where H.sub.0h(x,y)0 and:
h(x,y)H.sub.0.fwdarw.0 for {square root over (x.sup.2+y.sup.2)}.fwdarw.,
where H.sub.0 is a constant.

(43) It can be shown that the gravity field, g, of the infinitely extended domain can be represented by the following Cauchy-type integral:
g(r)=4.sub.0C.sup.(r,(z+H.sub.0)d.sub.z),
which can be re-written using matrix notation for the Cauchy integral from equation 2:

(44) $\begin{array}{cc}{g}_a\left(r^{}\right)=-{}_0{}_S{}_{}\left(z+H_0\right)\frac{r_{}-r_{}^{}}{{.Math.r-r^{}.Math.}^3}{n}_{}ds.& \left(3\right)\end{array}$

(45) Similarly, the gravity gradients can be expressed in matrix notation for the Cauchy integral from equation 2:

(46) $\begin{array}{cc}{g}_{\mathrm{av}}\left(r^{}\right)=\frac{g_a\left(r^{}\right)}{r_v^{}}=-{}_0{}_S{}_{}\frac{\left(z+H_0\right)}{{.Math.r-r^{}.Math.}^5}\left[3\left(r_v-r_v^{}\right)\left(r_{}-r_{}^{}\right)+{.Math.r-r^{}.Math.}^2{}_v\right]n_{}ds.& \left(4\right)\end{array}$

(47) We can provide explicit expressions for the components of the gravity and gravity gradients of the density contrast surface, taking into account the following relations for the components of the unit normal vector to the surface :

(48) $\begin{array}{cc}{n}_{x}ds=-\frac{h\left(x,y\right)}{x}dxdy=b_x\left(x,y\right)dxdy,n_{y}ds=-\frac{h\left(x,y\right)}{y}xy=b_y\left(x,y\right)dxdy,n_{z}ds=b_z\left(x,y\right)dxdy,\left(z+H_0\right)=h\left(x,y\right),& \left(5\right)\end{array}$
where:

(49) $b_x\left(x,y\right)=-\frac{h\left(x,y\right)}{x},b_y\left(x,y\right)=-\frac{h\left(x,y\right)}{y},b_z\left(x,y\right)=1.$
Substituting equations (5) into equations (3) and (4), we obtain:

(50) 0$\begin{array}{cc}{g}_a\left(r^{}\right)=-{}_0{}_S{}_{}h\left(x,y\right)\frac{{\tilde{r}}_{}-r_{}^{}}{{.Math.\tilde{r}-r^{}.Math.}^3}{b}_{}\left(x,y\right)dxdy,& \left(6\right)\end{array}$
for the gravity field, and:

(51) $\begin{array}{cc}{g}_{av}\left(r^{}\right)=-{}_0{}_S{}_{}\frac{h\left(x,y\right)}{{.Math.\tilde{r}-r^{}.Math.}^5}[3\left({\tilde{r}}_v-r_v^{}\right)\left({\tilde{r}}_{}-r_{}^{}\right)+{.Math.\tilde{r}-r^{}.Math.}^2{}_v]b_y\left(x,y\right)\mathrm{dx}\mathrm{dy}.& \left(7\right)\end{array}$
for the gravity gradients, where:
|{tilde over (r)}r|={square root over ((xx).sup.2+(yy).sup.2+(h(x,y)H.sub.0z).sup.2)},
{tilde over (r)}.sub.x=x,
{tilde over (r)}.sub.y=y,
{tilde over (r)}.sub.z=h(x,y)H.sub.0.

(52) For sedimentary basins, terrain (or bathymetry) corrections are generally calculated as the response due to the volume of earth bound by an upper surface of the digital elevation (or bathymetry) model and a lower surface of a plane at depth. To simulate sediment compaction and diagenesis causing a loss of porosity, densities are often parameterized from empirical observations using analytic density-depth functions obtained from with the terrain (or bathymetry) forming the datum. The variety of analytic density-depth functions in use span linear, quadratic, parabolic, exponential, hyperbolic, and polynomial functions. In general, the terrain can be estimated with a model having a vertically variable density:
=(z).

(53) It can be shown that the gravity field, g, of an infinitely extended domain can be represented as the following Cauchy-type integral:
g(r)=4.sub.0C.sup.(r[R(z)R(H.sub.0)]d.sub.z),(8)
where R(z) is the indefinite integral of density:
R(z)=(z)dz.
Equation (8) can be re-written using matrix notation for the Cauchy integral from equation 2:

(54) $\begin{array}{cc}{g}_{}\left(r^{}\right)=-{}_0{}_S{}_{}\left[R\left(z\right)-R\left(-H_0\right)\right]\frac{r_{}-r_{}^{}}{{.Math.r-r^{}.Math.}^3}{n}_{}ds.& \left(9\right)\end{array}$
Similarly, the gravity gradients can be expressed in matrix notation for the Cauchy integral from equation 2:

(55) $\begin{array}{cc}{g}_{\mathrm{av}}\left(r^{}\right)=\frac{g_a\left(r^{}\right)}{r_v^{}}={}_0{}_S{}_{}\frac{\left[R\left(z\right)-R\left(-H_0\right)\right]}{{.Math.r-r^{}.Math.}^5}\left[3\left(r_v-r_v^{}\right)\left(r_{}-r_{}^{}\right)+{.Math.r-r^{}.Math.}^2{}_v\right]n_{}ds.& \left(10\right)\end{array}$

(56) One can derive explicit expressions for the gravity fields and their gradients of the density contact surface, taking into account Equations (5) for the components of the unit normal vector to the surface, F:

(57) $\begin{array}{cc}{g}_{}\left(r^{}\right)=-{}_0{}_S{}_{}\left[R\left(z\right)-R\left(-H_0\right)\right]\frac{{\tilde{r}}_{}-r_{}^{}}{{.Math.\tilde{r}-r^{}.Math.}^3}{b}_{}ds,& \left(11\right)\end{array}$
and

(58) $\begin{array}{cc}{g}_v\left(r^{}\right)=\frac{g_{}\left(r^{}\right)}{r_v^{}}={}_0{}_S{}_{}\frac{\left[R\left(z\right)-R\left(-H_0\right)\right]}{{.Math.\tilde{r}-r^{}.Math.}^5}\left[3\left({\tilde{r}}_v-r_v^{}\right)\left({\tilde{r}}_{}-r_{}^{}\right)+{.Math.\tilde{r}-r^{}.Math.}^2{}_v\right]b_{}\left(x,y\right)ds.& \left(12\right)\end{array}$

(59) As an example of an embodiment of the method, one can consider a terrain with linear variations in the density with depth:
(z)=.sub.0+z,
such that:

(60) $R\left(z\right)={}_{0}z+\frac{1}{2}a{z}^2.$

(61) As another example of an embodiment of the method, one can consider a terrain with exponential variations in the density with depth:
(z)=.sub.0+exp(kz),
such that:

(62) $R\left(z\right)={}_{0}z+\frac{a}{k}\exp \left(\mathrm{kz}\right).$

(63) In one embodiment of this invention, the surface Cauchy-type integral equations 6 and 7 may be discretized by dividing the horizontal integration plane XY into a rectangular grid of N.sub.m cells with constant discretization of x and y in the x and y directions, respectively. In each cell P.sub.k (k=1, 2, . . . , N.sub.m), it is assumed that the corresponding part of the density contrast surface can be represented by an element of a plane described by the following equation:
z=h(x,y)H.sub.0=h.sup.(k)b.sub.x.sup.(k)(xx.sub.k)b.sub.y.sup.(k)(yy.sub.k)H.sub.0,
where (x,y) P.sub.k and (x.sub.k, y.sub.k) denotes the center of the cell P.sub.k. In this case:
b.sub.x(x,y)=b.sub.x.sup.(k),
b.sub.y(x,y)=b.sub.y.sup.(k),
b.sub.z(x,y)=b.sub.z.sup.(k)=1.
where (x,y) E P.sub.k. Equation 6 for the gravity field now may take the form:

(64) $\begin{array}{cc}{g}_{}\left(r^{}\right)=-{}_0\underset{k=1}{\overset{N_m}{.Math.}}{}_{P_k}{}_{}h\left(x,y\right)\frac{{\tilde{r}}_{}-r_{}^{}}{.Math.\tilde{r}-r^{}.Math.}{b}_{}^{\left(k\right)}dxdy.& \left(13\right)\end{array}$
Using the discrete model parameters introduced above, and discrete gravity data, g.sub.(r), we can represent the 3D forward modeling operator for the gravity field of the density contrast surface , equation 8, as:

(65) $\begin{array}{cc}{g}_a\left(r^{}\right)=\underset{k=1}{\overset{N_m}{.Math.}}{f}_{}^{\left(\mathrm{nk}\right)}{h}^{\left(k\right)}{b}_{}^{\left(k\right)},& \left(14\right)\end{array}$
where:

(66) 0$f_{}^{\left(\mathrm{nk}\right)}=-{}_0{}_{}\frac{{\tilde{r}}_{}^{\left(k\right)}-r_{}^{{\left(n\right)}^{}}}{{.Math.{\tilde{r}}_k-r_{}^{}.Math.}^3}xy,$
and:
|{tilde over (r)}r|={square root over ((x.sub.kx.sub.n).sup.2+(y.sub.ky.sub.n).sup.2+(h.sup.(k)H.sub.0z.sub.n).sup.2)},
{tilde over (r)}.sub.x.sup.(k)=x.sub.k,
{tilde over (r)}.sub.y.sup.(k)=y.sub.k,
{tilde over (r)}.sub.z.sup.(k)=h.sup.(k)H.sub.0.
{tilde over (r)}.sub.x.sup.(n)=x.sub.n,
{tilde over (r)}.sub.y.sup.(n)=y.sub.n,
{tilde over (r)}.sub.y.sup.(n)=y.sub.n,

(67) In a similar manner, the discrete form of equation 7 for the gravity gradients may be written as:

(68) $\begin{array}{cc}{g}_v\left(r^{}\right)={}_0\underset{k=1}{\overset{N_m}{.Math.}}{}_{P_k}{}_{}\frac{h\left(x,y\right)}{{.Math.\tilde{r}-r^{}.Math.}^5}[3\left({\tilde{r}}_v-r_v^{}\right)\left({\tilde{r}}_{}-r_{}^{}\right)+{.Math.\tilde{r}-r^{}.Math.}^2{}_v]b_{}^{\left(k\right)}dx\mathrm{dy}.& \left(15\right)\end{array}$
Using the discrete model parameters introduced above, and discrete gravity gradiometry data, g.sub.(r), we can represent the 3D forward modeling operator for the gravity gradients of the density contrast surface , equation 15, as:

(69) $\begin{array}{cc}{g}_v\left(r^{}\right)=\underset{k=1}{\overset{N_m}{.Math.}}{F}_v^{\left(\mathrm{nk}\right)}{h}^{\left(k\right)}{b}_{}^{\left(k\right)},& \left(16\right)\end{array}$
where:

(70) $F_v^{\left(\mathrm{nk}\right)}={}_0{}_{}\frac{1}{{.Math.{\tilde{r}}_k-r_{}^{}.Math.}^5}\left[3\left({\tilde{r}}_v-{\tilde{r}}_v^{\left(n\right)}\right)\left({\tilde{r}}_{}-{\tilde{r}}_{}^{\left(n\right)}\right)+{.Math.{\tilde{r}}_k-r_{}^{}.Math.}^2{}_v\right]xy,$

(71) In another embodiment of this invention, we may assume that the function h(x,y) describing the density contrast surface is given on a rectangular mesh :

(72) $.Math.=\left\{\left(x_i,y_j\right)\begin{array}{ccc}{x}_1=x^{}& x_I=x^{}& x_{i+1}=x_i+x\\ y_1=y^{}& y_J=y^{}& y_{j+1}=y_j+y\end{array}\right\}$
where i=1, 2, . . . , I and j=1, 2, . . . , J. In this case, the coordinates of the points located on the surface , {tilde over (r)}.sub.k={tilde over (r)}.sub.ij, are as follows:
{tilde over (r)}.sub.k=(x.sub.k,y.sub.k,h.sup.(k)H.sub.0)=(x.sub.i,y.sub.j,h.sup.(ij)H.sub.0)={tilde over (r)}.sub.ij,
where h.sup.(k)=h.sup.(ij)=h(x.sub.i,y.sub.j). Note that the nodes of the mesh are consecutively numbered along the horizontal and vertical directions. For the given numbering of the nodes, k=1, 2, . . . , K, where K=IJ, one can establish a simple one-to-one relationship between the index n and the double integer (i,j):
k=i+(j1)K.

(73) As illustrated in FIG. 7, we may divide each rectangular cell of the mesh into two triangles which would form the elementary cells P.sub.Lk (left triangular 44) and P.sub.Rk (right triangular 45). The corresponding equations of the plane parts, F.sub.Lk and F.sub.Rk, of the surface F, located just above the cells P.sub.Lk and P.sub.Rk can be written using the equations of the plane through the corresponding points {tilde over (r)}.sub.ij, {tilde over (r)}.sub.i+1,j, {tilde over (r)}.sub.i+1,j+1 and {tilde over (r)}.sub.ij, {tilde over (r)}.sub.i,j+1, {tilde over (r)}.sub.+i+1,j+1 for the left and right triangles, respectively:
(r{tilde over (r)}.sub.ij).Math.[({tilde over (r)}.sub.i+1,j{tilde over (r)}.sub.ij)({tilde over (r)}.sub.i+1,j+1{tilde over (r)}.sub.ij)]=0 (17)
for .sub.Lk, and:
(r{tilde over (r)}.sub.ij).Math.[({tilde over (r)}.sub.i,j+1{tilde over (r)}.sub.ij)({tilde over (r)}.sub.i+1,j+1{tilde over (r)}.sub.ij)]=0 (18)
for .sub.Rk. One can write an equivalent form of equation (12) as:
.sub.Lx.sup.(k)x+.sub.Ly.sup.(k)y+.sub.Lz.sup.(k)z+.sub.L0.sup.(k)=0,
where:

(74) $a_{\mathrm{Lx}}^{\left(k\right)}=.Math.\begin{array}{ccc}1& y_{\mathrm{ij}}& z_{\mathrm{ij}}\\ 1& y_{i+1,j}& z_{i+1,j}\\ 1& y_{i+1,j+1}& z_{i+1,j+1}\end{array}.Math.$$a_{\mathrm{Ly}}^{\left(k\right)}=.Math.\begin{array}{ccc}{x}_{\mathrm{ij}}& 1& z_{\mathrm{ij}}\\ x_{i+1,j}& 1& z_{i+1,j}\\ x_{i+1,j+1}& 1& z_{i+1,j+1}\end{array}.Math.$$a_{\mathrm{Lz}}^{\left(k\right)}=.Math.\begin{array}{ccc}{x}_{\mathrm{ij}}& y_{\mathrm{ij}}& 1\\ x_{i+1,j}& y_{i+1,j}& 1\\ x_{i+1,j+1}& y_{i+1,j+1}& 1\end{array}.Math.$$a_{L0}^{\left(k\right)}=.Math.\begin{array}{ccc}{x}_{\mathrm{ij}}& y_{\mathrm{ij}}& z_{\mathrm{ij}}\\ x_{i+1,j}& y_{i+1,j}& z_{i+1,j}\\ x_{i+1,j+1}& y_{i+1,j+1}& z_{i+1,j+1}\end{array}.Math.$
Similarly, one may write an equivalent form of equation (18) as:
.sub.Rx.sup.(k)x+.sub.Ry.sup.(k)y+.sub.Rz.sup.(k)z+.sub.R0.sup.(k)=0,
where:

(75) $a_{\mathrm{Rx}}^{\left(k\right)}=.Math.\begin{array}{ccc}1& y_{\mathrm{ij}}& z_{\mathrm{ij}}\\ 1& y_{i,j+1}& z_{i,j+1}\\ 1& y_{i+1,j+1}& z_{i+1,j+1}\end{array}.Math.$$a_{\mathrm{Ry}}^{\left(k\right)}=.Math.\begin{array}{ccc}{x}_{\mathrm{ij}}& 1& z_{\mathrm{ij}}\\ x_{i,j+1}& 1& z_{i,j+1}\\ x_{i+1,j+1}& 1& z_{i+1,j+1}\end{array}.Math.$$a_{\mathrm{Rz}}^{\left(k\right)}=.Math.\begin{array}{ccc}{x}_{\mathrm{ij}}& y_{\mathrm{ij}}& 1\\ x_{i,j+1}& y_{i,j+1}& 1\\ x_{i+1,j+1}& y_{i+1,j+1}& 1\end{array}.Math.$$a_{R0}^{\left(k\right)}=.Math.\begin{array}{ccc}{x}_{\mathrm{ij}}& y_{\mathrm{ij}}& z_{\mathrm{ij}}\\ x_{i,j+1}& y_{i,j+1}& z_{i,j+1}\\ x_{i+1,j+1}& y_{i+1,j+1}& z_{i+1,j+1}\end{array}.Math.$
One may introduce the following notations:
h.sub.L.sup.(k)=[.sub.Lx.sup.(k)x.sub.Lk+.sub.Ly.sup.(k)y.sub.Lk+.sub.L0.sup.(k)]/.sub.Lz.sup.(k), (19)
h.sub.R.sup.(k)=[.sub.Rx.sup.(k)x.sub.Rk+.sub.Ry.sup.(k)y.sub.Rk+.sub.R0.sup.(k)]/.sub.Rz.sup.(k), (20)
Where (x.sub.Lk, y.sub.Lk) and (x.sub.Rk,y.sub.Rk) denote the center of the cells P.sub.Lk and P.sub.Rk respectively. Substituting equations (19) and (20) into equation (5), one obtains:
b.sub.Lx.sup.(k)=.sub.Lx.sup.(k)/.sub.Lz.sup.(k),
b.sub.Ly.sup.(k)=.sub.Ly.sup.(k)/.sub.Lz.sup.(k),
b.sub.Rx.sup.(k)=.sub.Rx.sup.(k)/.sub.Rz.sup.(k),
b.sub.Ry.sup.(k)=.sub.Ry.sup.(k)/.sub.Rz.sup.(k).

(76) Using the discrete model parameters introduced above, and discrete gravity data, g.sub.(r), we may represent the 3D forward modeling operator for the gravity field of the density contrast surface , equation 8, as:

(77) $\begin{array}{cc}{g}_{}\left(r^{}\right)=\underset{k=1}{\overset{N_m}{.Math.}}{f}_L^{\left(\mathrm{nk}\right)}{h}_L^{\left(k\right)}{b}_L^{\left(k\right)}+\underset{k=1}{\overset{N_m}{.Math.}}{f}_R^{\left(\mathrm{nk}\right)}{h}_R^{\left(k\right)}{b}_R^{\left(k\right)},& \left(21\right)\end{array}$
where

(78) $f_L^{\left(\mathrm{nk}\right)}=-{}_0{}_{}\frac{{\tilde{r}}_L^{\left(k\right)}-r_{}^{\left(n\right)}}{{.Math.{\tilde{r}}_{\mathrm{Lk}}-r_{}^{}.Math.}^3}xy,f_R^{\left(\mathrm{nk}\right)}=-{}_0{}_{}\frac{{\tilde{r}}_R^{\left(k\right)}-r_{}^{\left(n\right)}}{{.Math.{\tilde{r}}_{\mathrm{Rk}}-r_{}^{}.Math.}^3}xy,$
and:
{tilde over (r)}.sub.Lx.sup.(k)=x.sub.Lk,
{tilde over (r)}.sub.Ly.sup.(k)=y.sub.Lk,
{tilde over (r)}.sub.Lz.sup.(k)=h.sub.L.sup.(k)H.sub.0,
{tilde over (r)}.sub.Rx.sup.(k)=x.sub.Rk,
{tilde over (r)}.sub.Ry.sup.(k)=y.sub.Rk,
{tilde over (r)}.sub.Rz.sup.(k)=h.sub.R.sup.(k)=H.sub.0,
{tilde over (r)}.sub.x.sup.(n)=x.sub.n,
{tilde over (r)}.sub.y.sup.(n)=y.sub.n,
{tilde over (r)}.sub.z.sup.(n)=z.sub.n.

(79) In a similar manner, one may derive expressions for the gravity gradients using equation (16).

(80) The Cauchy integrals are evaluated independently for each data at each observation location. In one embodiment of the method, the Cauchy integrals can be evaluated on a grid with increasing discretization that is a function of distance from the observation location, and which may terminate at some distance from the observation location that is determined from the integrated sensitivity of the data with respect to the 3D density model.

EXAMPLE 1

(81) FIG. 8 is a plan view of five anomalous density bodies, including: a tabular block 46 of 0.20 g/cm.sup.3 density contrast buried at 500 m depth with a strike of 150, strike length of 600 m, width of 200 m, and depth extent of 400 m; a dipping sheet 47 of 0.15 g/cm.sup.3 density contrast buried at 0 m depth with a dip of 75 m, dip of 45, strike of 50, strike length of 3000 m, and down-dip depth extent of 1000 m; a sphere 48 of 0.30 g/cm.sup.3 density contrast buried at 50 m depth with a radius of 150 m; a truncated cone 49 that narrows towards its base of 0.40 g/cm.sup.3 density contrast buried at 150 m depth with an upper radius of 100 m, a lower radius of 70 m, and a depth extent of 200 m; and a tunnel 50 of density contrast 2.7 g/cm.sup.3 buried at 3 m depth, with a width of 4 m, depth extent of 4 m, and a length of 500 m.

(82) FIG. 9 is a plan view of the digital elevation model as measured from a LiDAR survey over the synthetic survey area.

(83) FIG. 10 is a plan view of the digital elevation model as produced from the merging of the LiDAR digital elevation model over the synthetic survey area and the SRTM digital elevation model over a 10 km area beyond the synthetic survey area.

(84) FIG. 11 is the free-air vertical gravity gradient (Gzz) data simulated with 1 E/Hz noise added for the synthetic density model and digital elevation model with a homogeneous terrain density of 2.67 g/cm.sup.3.

(85) FIG. 12 is the terrain corrected vertical gravity gradient (Gzz) data obtained from terrain correction of the synthetic free air vertical gravity gradient (Gzz) data and merged LiDAR and SRTM digital elevation models with a homogeneous terrain density of 2.67 g/cm.sup.3.

(86) Embodiments of the present invention may comprise or utilize a special purpose or general-purpose computer including computer hardware, as discussed in greater detail below. Embodiments within the scope of the present invention also include physical and other computer-readable media for carrying or storing computer-executable instructions and/or data structures. Such computer-readable media can be any available media that can be accessed by a general purpose or special purpose computer system. Computer-readable media that store computer-executable instructions are physical non-transitory storage media. Computer-readable media that carry computer-executable instructions are transmission media. Thus, by way of example, and not limitation, embodiments of the invention can comprise at least two distinctly different kinds of computer-readable media: physical non-transitory storage media and transmission media.

(87) Physical non-transitory storage media includes RAM, ROM, EEPROM, CD-ROM or other optical disk storage, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store desired program code means in the form of computer-executable instructions or data structures and which can be accessed by a general purpose or special purpose computer.

(88) A network is defined as one or more data links that enable the transport of electronic data between computer systems and/or modules and/or other electronic devices. When information is transferred or provided over a network or another communications connection (either hardwired, wireless, or a combination of hardwired or wireless) to a computer, the computer properly views the connection as a transmission medium. Transmissions media can include a network and/or data links which can be used to carry or desired program code means in the form of computer-executable instructions or data structures and which can be accessed by a general purpose or special purpose computer. Combinations of the above should also be included within the scope of computer-readable media.

(89) Further, upon reaching various computer system components, program code means in the form of computer-executable instructions or data structures can be transferred automatically from transmission media to physical storage media (or vice versa). For example, computer-executable instructions or data structures received over a network or data link can be buffered in RAM within a network interface module (e.g., a NIC), and then eventually transferred to computer system RAM and/or to less volatile physical storage media at a computer system. Thus, it should be understood that physical storage media can be included in computer system components that also (or even primarily) utilize transmission media.

(90) Computer-executable instructions comprise, for example, instructions and data which cause a general purpose computer, special purpose computer, or special purpose processing device to perform a certain function or group of functions. The computer executable instructions may be, for example, binaries, intermediate format instructions such as assembly language, or even source code. Although the subject matter has been described in language specific to structural features and/or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the described features or acts described above. Rather, the described features and acts are disclosed as example forms of implementing the claims.

(91) Those skilled in the art will appreciate that the invention may be practiced in network computing environments with many types of computer system configurations, including, personal computers, desktop computers, laptop computers, message processors, hand-held devices, multi-processor systems, microprocessor-based or programmable consumer electronics, network PCs, minicomputers, mainframe computers, mobile telephones, PDAs, pagers, routers, switches, and the like. The invention may also be practiced in distributed system environments where local and remote computer systems, which are linked (either by hardwired data links, wireless data links, or by a combination of hardwired and wireless data links) through a network, both perform tasks. In a distributed system environment, program modules may be located in both local and remote memory storage devices.

(92) The methods disclosed herein comprise one or more steps or actions for achieving the described method. The method steps and/or actions may be interchanged with one another without departing from the scope of the present invention. In other words, unless a specific order of steps or actions is required for proper operation of the embodiment, the order and/or use of specific steps and/or actions may be modified without departing from the scope of the present invention.

(93) While specific embodiments and applications of the present invention have been illustrated and described, it is to be understood that the invention is not limited to the precise configuration and components disclosed herein. Various modifications, changes, and variations which will be apparent to those skilled in the art may be made in the arrangement, operation, and details of the methods and systems of the present invention disclosed herein without departing from the spirit and scope of the invention.