Method of terrain correction for potential field geophysical survey data

[0001] This application claims the benefit of and priority to U.S. Provisional Application No. 61/578,571, filed Dec. 21, 2011, which is incorporated herein by reference in its entirety.

[0002] This application hereby incorporates the following publication by reference in its entirety: Zhdanov, M. S., 1988, Integral transforms in geophysics: Springer-Verlag, Berlin.

[0003] 1. The Field of the Invention

[0004] The present disclosure relates in general to methods of terrain correction for potential field geophysical survey data, particularly gravity and gravity gradiometry data.

[0005] 2. The Related Technology

[0006] Potential field geophysical surveys are performed by measuring vector and/or tensors of a potential field, i.e., the first or second order spatial derivatives of the magnetic and/or gravity potentials. Potential field geophysical surveys are widely used for mapping subsurface geological and man-made structures in mineral, hydrocarbon, geothermal and groundwater exploration, and hydrocarbon, geothermal and groundwater resource monitoring, earth systems modeling, and tunnel and underground facility (UGF) detection.

[0007] Potential field geophysical data usually include at least one component of the gravity and/or magnetic vector and/or tensor fields. Measured potential field data contains the linear superposition of responses from both the terrain, and subsurface geological structures located beneath it.

[0008] To provide economical reconnaissance of subsurface geological structures, potential field sensors are often deployed from moving platforms such as airplanes, helicopters, airships, unattended aerial vehicles (UAV), borehole logging instruments, vessels, submarines, and vehicles.

[0009] For example, airborne potential field surveys from fixed wing aircraft typically acquire 500 line km of data each day, and airborne potential field surveys from helicopters typically acquire 200 line km of data each day. Airborne potential field surveys typically contain multiple survey lines draped over a variable terrain surface that aggregate as hundreds to thousands of line kilometers of potential field data measured every few meters and cover an area hundreds to thousands of square kilometers in size. The airborne platform's drape surface is constrained by the minimum ground clearance that the airborne platform is permitted to fly, and the maximum rate of climb or descent of the airborne platform relative to the terrain.

[0010] In practice, terrain measurements are usually embodied in Digital Elevation Models (DEM) as Easting (or Longitude), Northing (or Latitude) and Elevation constructed from the processing of (Differential) Global Position Systems ((D)GPS) and/or Radio Detection and Ranging (RADAR) and/or Laser Imaging Detection and Ranging (LIDAR) and/or Synthetic Aperture Radar (SAR) data measured from airborne platforms, or from pubic geospatial databases such as NASA's Shuttle Topography Radar Mission (SRTM).

[0011] The effect of terrain contributes to the measured potential field data. The effect of terrain actually dominates free-air gravity and gravity gradiometry data. To give an indication of the magnitude of the different responses in gravity gradiometry, a typical terrain response maybe in the order of hundreds of Eotvos (Eö), whereas the target response maybe in the order of tens of Eö.

[0012] After the potential field geophysical data have been acquired, a terrain correction needs to applied to compensate for the variations in topography (White et al., 1995, U.S. Pat. No. 5,402,340).

[0013] One method for terrain correction involves three-dimensional (3D) modeling of the potential field response in the spatial domain due to the DEM discretized using a regular or irregular grid of right rectangular prisms (“cells”) bound above by the DEM, and below by a plane that passes through the lowest elevation of the DEM. Expressions for the potential field response of each cell contain volume integrals, and these may be evaluated using analytical or numerical techniques. The cells in the grid are usually populated with a constant physical property, e.g., a constant density of 2.67 g/cc corresponding to an average upper continental crust density. In some implementations, for example for sedimentary basins, the grid is populated with a spatially variable physical property, e.g., a density analytically increasing with depth corresponding to sediment compaction. To reduce computational requirements, the cell size is increased with increasing distance from the observation point (Brewster, 2006, US Patent Application US2006/0036367 A1; Davies, 2009, US Patent Application US2009/0252372 A1; Barnes, 2009, US Patent Application US2009/0287464 A1; Barnes et al., 2010, US Patent Application US 2010/0094556 A1; Davies, 2010, US Patent Application US 2010/0206557).

[0014] Another method for terrain correction involves three-dimensional (3D) modeling of the potential field response in the Fourier (or wavenumber) domain due to the DEM discretized using a polygon bound above by the DEM, and below by a plane that passes through the lowest elevation of the DEM. Expressions for the potential field response of each cell contain Fourier transforms, and these may be evaluated using numerical techniques such as Fast Fourier Transforms (FFTs). The polygon is usually populated with a constant physical property, e.g., a constant density of 2.67 g/cc corresponding to an average upper continental crust density. In some implementations, for example for sedimentary basins, the polygon is populated with a spatially variable physical property, e.g., a density analytically increasing with depth corresponding to sediment compaction.

[0015] The accuracies of these terrain corrections are limited by numerical considerations, such as the discretization of the cells, accuracy of the volume integration, and accuracy of the Fourier transforms used. Accurate terrain corrections are computationally time consuming.

[0016] Despite improved navigation, instrumentation, and data processing techniques, there remains a need for improved terrain correction techniques for potential field geophysical data. This will be more important as next generation of potential field instruments such as very low noise (e.g., 1 Eö/√Hz) gravity gradiometers are currently being developed and tested.

[0017] The embodiments disclosed herein are related to systems, methods, and computer readable medium for terrain correction for potential field geophysical data measured from at least one potential field sensor at at least one measurement position based on representation of the terrain response as a 3D Cauchy-type surface integral over the topography of the earth surface. In the systems, methods, and computer readable medium at least one component of gravity and/or magnetic scalar, vector, and/or tensor data with at least one potential field sensor in at least one measurement position along at least one survey line is measured. The terrain as a digital elevation model is also measured. An appropriate physical property function to describe the physical property distribution within the terrain is selected. The terrain response from the digital elevation model with a physical property function for the at least one component of gravity and/or magnetic vector and/or tensor data with at least one potential field sensor in at least one measurement position along at least one survey line, using 3D Cauchy-type surface integral over the terrain surface is forward modeled. The system terrain response for the at least one component of gravity and/or magnetic vector and/or tensor data with at least one potential field sensor in at least one measurement position along at least one survey line is calculated. The terrain corrected potential field data is also calculated.

[0031] 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.

[0032] 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).

[0033] Measured potential field data contains the linear superposition of responses from both the terrain and the subsurface geological structures located beneath it.

[0034] 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.

[0035] 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.

[0036] 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.

[0037] 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.

[0038] 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.

[0039] 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.

[0040] 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.

[0041] 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.

[0042] At least one embodiment of this method can be used in geophysical exploration for mineral, hydrocarbon, geothermal, and groundwater resources.

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

[0044] At least one embodiment of this method can be used in geophysical monitoring for earth systems.

[0045] At least one embodiment of this method can be used for tunnels and UGF detection.

[0046] 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.

[0047] 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.

[0048] 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.

[0049] 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.

[0050] 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.

[0051] 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.

[0052] The concept of a three-dimensional Cauchy-type integral was introduced by Zhdanov, 1988, and is represented by the following expression:

[0053] $\begin{array}{cc}{C}^S\left(r,\phi \right)=-\frac{1}{4\pi }\int {\int }_S\left[\left(n·\phi \right)\nabla \frac{1}{r-r^{\prime }}+\left(n\times \phi \right)\times \nabla \frac{1}{r-r^{\prime }}\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^S(r, φ).

[0054] The Cauchy-type integral C^S(r, φ) can be represented in matrix notation. The vectors C^S, φ, n, and

[0055] $\nabla \frac{1}{r-r^{\prime }}$
are presented in some Cartesian basis {d_x, d_y, d_z}, where the z axis is directed upward, as:

[0056] $C^S=C_{\alpha }^S{d}_{\alpha }$ $\phi ={\phi }_{\beta }{d}_{\beta }$ $n=n_{\gamma }{d}_{\gamma }$ $\nabla \frac{1}{r-r^{\prime }}=-\frac{r-r^{\prime }}{{r-r^{\prime }}^3}=-\frac{r_{\eta }-r_{\eta }^{\prime }}{{r-r^{\prime }}^3}{d}_{\eta }$
where r_η=η 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:

[0057] $\begin{array}{cc}{C}^S\left(r,\phi \right)=-\frac{1}{4\pi }\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}{\phi }_{\beta }\frac{r_{\eta }-r_{\eta }^{\prime }}{{r-r^{\prime }}^3}{n}_{\gamma }ds,& \left(2\right)\end{array}$
where the four-index Δ-symbol is expressed in terms of the symmetric Kronecker symbol δ_αβ:

[0058] ${\delta }_{\mathrm{\alpha \beta }}={\delta }_{\mathrm{\alpha \beta }}{\delta }_{\mathrm{\gamma \eta }}+{\delta }_{\mathrm{\alpha \eta }}{\delta }_{\mathrm{\beta \gamma }}-{\delta }_{\mathrm{\alpha \gamma }}{\delta }_{\mathrm{\beta \eta }},{\delta }_{\mathrm{\alpha \beta }}=\{\begin{array}{cc}1,& \alpha =\beta ,\\ 0,& \alpha \ne \beta .\end{array}$

[0059] 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₀, and a horizontal plane P 43, z=h₀, where H₀≥h(x,y)≥0 and:
h(x,y)−H₀→0 for √{square root over (x²+y²)}→∞,
where H₀ is a constant.

[0060] 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πγρ₀C^Γ(r′,(z+H₀)d_z),
which can be re-written using matrix notation for the Cauchy integral from equation 2:

[0061] $\begin{array}{cc}{g}_a\left(r^{\prime }\right)=-{\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\left(z+H_0\right)\frac{r_{\eta }-r_{\eta }^{\prime }}{{r-r^{\prime }}^3}{n}_{\gamma }ds.& \left(3\right)\end{array}$

[0062] Similarly, the gravity gradients can be expressed in matrix notation for the Cauchy integral from equation 2:

[0063] $\begin{array}{cc}{g}_{\mathrm{av}}\left(r^{\prime }\right)=\frac{\partial g_a\left(r^{\prime }\right)}{\partial r_v^{\prime }}=-{\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{\left(z+H_0\right)}{{r-r^{\prime }}^5}\left[3\left(r_v-r_v^{\prime }\right)\left(r_{\eta }-r_{\eta }^{\prime }\right)+{r-r^{\prime }}^2{\delta }_{v\eta }\right]n_{\gamma }ds.& \left(4\right)\end{array}$

[0064] 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 Γ:

[0065] $\begin{array}{cc}{n}_{x}ds=-\frac{\partial h\left(x,y\right)}{\partial x}dxdy=b_x\left(x,y\right)dxdy,n_{y}ds=-\frac{\partial h\left(x,y\right)}{\partial y}dxdy=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:

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

[0067] $\begin{array}{cc}{g}_a\left(r^{\prime }\right)=-{\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}h\left(x,y\right)\frac{{\tilde{r}}_{\eta }-r_{\eta }^{\prime }}{{\tilde{r}-r^{\prime }}^3}{b}_{\gamma }\left(x,y\right)dxdy,& \left(6\right)\end{array}$
for the gravity field, and:

[0068] $\begin{array}{cc}{g}_{av}\left(r^{\prime }\right)=-{\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{h\left(x,y\right)}{{\tilde{r}-r^{\prime }}^5}[3\left({\tilde{r}}_v-r_v^{\prime }\right)\left({\tilde{r}}_{\eta }-r_{\eta }^{\prime }\right)+{\tilde{r}-r^{\prime }}^2{\delta }_{v\eta }]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 ((x−x′)²+(y−y′)²+(h(x,y)−H₀−z′)²)},
{tilde over (r)}_x=x,
{tilde over (r)}_y=y,
{tilde over (r)}_z=h(x,y)−H₀.

[0069] 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).

[0070] 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πγρ₀C^Γ(r′[R(z)−R(−H₀)]d_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:

[0071] $\begin{array}{cc}{g}_{\alpha }\left(r^{\prime }\right)=-{\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\left[R\left(z\right)-R\left(-H_0\right)\right]\frac{r_{\eta }-r_{\eta }^{\prime }}{{r-r^{\prime }}^3}{n}_{\gamma }ds.& \left(9\right)\end{array}$
Similarly, the gravity gradients can be expressed in matrix notation for the Cauchy integral from equation 2:

[0072] $\begin{array}{cc}{g}_{\mathrm{av}}\left(r^{\prime }\right)=\frac{\partial g_a\left(r^{\prime }\right)}{\partial r_v^{\prime }}={\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{\left[R\left(z\right)-R\left(-H_0\right)\right]}{{r-r^{\prime }}^5}\left[3\left(r_v-r_v^{\prime }\right)\left(r_{\eta }-r_{\eta }^{\prime }\right)+{r-r^{\prime }}^2{\delta }_{v\eta }\right]n_{\gamma }ds.& \left(10\right)\end{array}$

[0073] 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:

[0074] $\begin{array}{cc}{g}_{\alpha }\left(r^{\prime }\right)=-{\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\left[R\left(z\right)-R\left(-H_0\right)\right]\frac{{\tilde{r}}_{\eta }-r_{\eta }^{\prime }}{{\tilde{r}-r^{\prime }}^3}{b}_{\gamma }ds,& \left(11\right)\end{array}$
and

[0075] $\begin{array}{cc}{g}_{\alpha v}\left(r^{\prime }\right)=\frac{\partial g_{\alpha }\left(r^{\prime }\right)}{\partial r_v^{\prime }}={\mathrm{\gamma \rho }}_0\int {\int }_S{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{\left[R\left(z\right)-R\left(-H_0\right)\right]}{{\tilde{r}-r^{\prime }}^5}\left[3\left({\tilde{r}}_v-r_v^{\prime }\right)\left({\tilde{r}}_{\eta }-r_{\eta }^{\prime }\right)+{\tilde{r}-r^{\prime }}^2{\delta }_{v\eta }\right]b_{\gamma }\left(x,y\right)ds.& \left(12\right)\end{array}$

[0076] As an example of an embodiment of the method, one can consider a terrain with linear variations in the density with depth:
ρ(z)=ρ₀+αz,
such that:

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

[0078] As another example of an embodiment of the method, one can consider a terrain with exponential variations in the density with depth:
ρ(z)=ρ₀+αexp(kz),
such that:

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

[0080] 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_m cells with constant discretization of Δx and Δy in the x and y directions, respectively. In each cell P_k (k=1, 2, . . . , N_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₀=h^(k)−b_x^(k)(x−x_k)−b_y^(k)(y−y_k)−H₀,
where (x,y) ϵ P_k and (x_k, y_k) denotes the center of the cell P_k. In this case:
b_x(x,y)=b_x^(k),
b_y(x,y)=b_y^(k),
b_z(x,y)=b_z^(k)=1.
where (x,y) E P_k. Equation 6 for the gravity field now may take the form:

[0081] $\begin{array}{cc}{g}_{\alpha }\left(r^{\prime }\right)=-{\mathrm{\gamma \rho }}_0\sum _{k=1}^{N_m}{\int }_{P_k}\int {\Delta }_{\mathrm{\alpha \beta \gamma \eta }}h\left(x,y\right)\frac{{\tilde{r}}_{\eta }-r_{\eta }^{\prime }}{\tilde{r}-r^{\prime }}{b}_{\gamma }^{\left(k\right)}dxdy.& \left(13\right)\end{array}$
Using the discrete model parameters introduced above, and discrete gravity data, g_α(r′), we can represent the 3D forward modeling operator for the gravity field of the density contrast surface Γ, equation 8, as:

[0082] $\begin{array}{cc}{g}_a\left(r^{\prime }\right)=\sum _{k=1}^{N_m}{f}_{\mathrm{\alpha \gamma }}^{\left(\mathrm{nk}\right)}{h}^{\left(k\right)}{b}_{\gamma }^{\left(k\right)},& \left(14\right)\end{array}$
where:

[0083] $f_{\mathrm{\alpha \gamma }}^{\left(\mathrm{nk}\right)}=-{\mathrm{\gamma \rho }}_0{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{{\tilde{r}}_{\eta }^{\left(k\right)}-r_{\eta }^{{\left(n\right)}^{\prime }}}{{{\tilde{r}}_k-r_{\eta }^{\prime }}^3}\Delta x\Delta y,$
and:
|{tilde over (r)}−r′|=√{square root over ((x_k−x_n′)²+(y_k−y_n′)²+(h^(k)−H₀−z_n′)²)},
{tilde over (r)}_x^(k)=x_k,
{tilde over (r)}_y^(k)=y_k,
{tilde over (r)}_z^(k)=h^(k)−H₀.
{tilde over (r)}_x⁽ⁿ⁾′=x_n′,
{tilde over (r)}_y⁽ⁿ⁾′=y_n′,
{tilde over (r)}_y⁽ⁿ⁾′=y_n′,

[0084] In a similar manner, the discrete form of equation 7 for the gravity gradients may be written as:

[0085] $\begin{array}{cc}{g}_{\alpha v}\left(r^{\prime }\right)={\mathrm{\gamma \rho }}_0\sum _{k=1}^{N_m}{\int }_{P_k}\int {\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{h\left(x,y\right)}{{\tilde{r}-r^{\prime }}^5}[3\left({\tilde{r}}_v-r_v^{\prime }\right)\left({\tilde{r}}_{\eta }-r_{\eta }^{\prime }\right)+{\tilde{r}-r^{\prime }}^2{\delta }_{v\eta }]b_{\gamma }^{\left(k\right)}dx\mathrm{dy}.& \left(15\right)\end{array}$
Using the discrete model parameters introduced above, and discrete gravity gradiometry data, g_αν(r′), we can represent the 3D forward modeling operator for the gravity gradients of the density contrast surface Γ, equation 15, as:

[0086] $\begin{array}{cc}{g}_{\alpha v}\left(r^{\prime }\right)=\sum _{k=1}^{N_m}{F}_{\alpha v\gamma }^{\left(\mathrm{nk}\right)}{h}^{\left(k\right)}{b}_{\gamma }^{\left(k\right)},& \left(16\right)\end{array}$
where:

[0087] $F_{\alpha v\gamma }^{\left(\mathrm{nk}\right)}={\mathrm{\gamma \rho }}_0{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{1}{{{\tilde{r}}_k-r_{\eta }^{\prime }}^5}\left[3\left({\tilde{r}}_v-{\tilde{r}}_v^{\left(n\right)\prime }\right)\left({\tilde{r}}_{\eta }-{\tilde{r}}_{\eta }^{\left(n\right)\prime }\right)+{{\tilde{r}}_k-r_{\eta }^{\prime }}^2{\delta }_{v\eta }\right]\Delta x\Delta y,$

[0088] 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 Σ:

[0089] $\sum =\left\{\left(x_i,y_j\right)\begin{array}{ccc}{x}_1=x^{\prime }& x_I=x^{″}& x_{i+1}=x_i+\Delta x\\ y_1=y^{\prime }& y_J=y^{″}& y_{j+1}=y_j+\Delta 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)}_k={tilde over (r)}_ijϵΓ, are as follows:
{tilde over (r)}_k=(x_k,y_k,h^(k)−H₀)=(x_i,y_j,h^(ij)−H₀)={tilde over (r)}_ij,
where h^(k)=h^(ij)=h(x_i,y_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+(j−1)K.

[0090] As illustrated in FIG. 7, we may divide each rectangular cell of the mesh into two triangles which would form the elementary cells P_Lk (left triangular 44) and P_Rk (right triangular 45). The corresponding equations of the plane parts, F_Lk and F_Rk, of the surface F, located just above the cells P_Lk and P_Rk can be written using the equations of the plane through the corresponding points {tilde over (r)}_ij, {tilde over (r)}_i+1,j, {tilde over (r)}_i+1,j+1 and {tilde over (r)}_ij, {tilde over (r)}_i,j+1, {tilde over (r)}_+i+1,j+1 for the left and right triangles, respectively:
(r−{tilde over (r)}_ij)·[({tilde over (r)}_i+1,j−{tilde over (r)}_ij)×({tilde over (r)}_i+1,j+1−{tilde over (r)}_ij)]=0   (17)
for Γ_Lk, and:
(r−{tilde over (r)}_ij)·[({tilde over (r)}_i,j+1−{tilde over (r)}_ij)×({tilde over (r)}_i+1,j+1−{tilde over (r)}_ij)]=0   (18)
for Γ_Rk. One can write an equivalent form of equation (12) as:
α_Lx^(k)x+α_Ly^(k)y+α_Lz^(k)z+α_L0^(k)=0,
where:

[0091] $a_{\mathrm{Lx}}^{\left(k\right)}=\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}$ $a_{\mathrm{Ly}}^{\left(k\right)}=\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}$ $a_{\mathrm{Lz}}^{\left(k\right)}=\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}$ $a_{L0}^{\left(k\right)}=\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}$
Similarly, one may write an equivalent form of equation (18) as:
α_Rx^(k)x+α_Ry^(k)y+α_Rz^(k)z+α_R0^(k)=0,
where:

[0092] $a_{\mathrm{Rx}}^{\left(k\right)}=\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}$ $a_{\mathrm{Ry}}^{\left(k\right)}=\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}$ $a_{\mathrm{Rz}}^{\left(k\right)}=\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}$ $a_{R0}^{\left(k\right)}=\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}$
One may introduce the following notations:
h_L^(k)=−[α_Lx^(k)x_Lk+α_Ly^(k)y_Lk+α_L0^(k)]/α_Lz^(k),   (19)
h_R^(k)=−[α_Rx^(k)x_Rk+α_Ry^(k)y_Rk+α_R0^(k)]/α_Rz^(k),   (20)
Where (x_Lk, y_Lk) and (x_Rk,y_Rk) denote the center of the cells P_Lk and P_Rk respectively. Substituting equations (19) and (20) into equation (5), one obtains:
b_Lx^(k)=α_Lx^(k)/α_Lz^(k),
b_Ly^(k)=α_Ly^(k)/α_Lz^(k),
b_Rx^(k)=α_Rx^(k)/α_Rz^(k),
b_Ry^(k)=α_Ry^(k)/α_Rz^(k).

[0093] Using the discrete model parameters introduced above, and discrete gravity data, g_α(r′), we may represent the 3D forward modeling operator for the gravity field of the density contrast surface Γ, equation 8, as:

[0094] $\begin{array}{cc}{g}_{\alpha }\left(r^{\prime }\right)=\sum _{k=1}^{N_m}{f}_{L\mathrm{\alpha \gamma }}^{\left(\mathrm{nk}\right)}{h}_L^{\left(k\right)}{b}_{L\gamma }^{\left(k\right)}+\sum _{k=1}^{N_m}{f}_{R\mathrm{\alpha \gamma }}^{\left(\mathrm{nk}\right)}{h}_R^{\left(k\right)}{b}_{R\gamma }^{\left(k\right)},& \left(21\right)\end{array}$
where

[0095] $f_{L\mathrm{\alpha \gamma }}^{\left(\mathrm{nk}\right)}=-{\mathrm{\gamma \rho }}_0{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{{\tilde{r}}_{L\eta }^{\left(k\right)}-r_{\eta }^{\left(n\right)\prime }}{{{\tilde{r}}_{\mathrm{Lk}}-r_{\eta }^{\prime }}^3}\Delta x\Delta y,f_{R\alpha \gamma }^{\left(\mathrm{nk}\right)}=-{\mathrm{\gamma \rho }}_0{\Delta }_{\mathrm{\alpha \beta \gamma \eta }}\frac{{\tilde{r}}_{R\eta }^{\left(k\right)}-r_{\eta }^{\left(n\right)\prime }}{{{\tilde{r}}_{\mathrm{Rk}}-r_{\eta }^{\prime }}^3}\Delta x\Delta y,$
and:
{tilde over (r)}_Lx^(k)=x_Lk,
{tilde over (r)}_Ly^(k)=y_Lk,
{tilde over (r)}_Lz^(k)=h_L^(k)−H₀,
{tilde over (r)}_Rx^(k)=x_Rk,
{tilde over (r)}_Ry^(k)=y_Rk,
{tilde over (r)}_Rz^(k)=h_R^(k)=H₀,
{tilde over (r)}_x⁽ⁿ⁾′=x_n′,
{tilde over (r)}_y⁽ⁿ⁾′=y_n′,
{tilde over (r)}_z⁽ⁿ⁾′=z_n′.

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

[0097] 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.

[0098] FIG. 8 is a plan view of five anomalous density bodies, including: a tabular block 46 of 0.20 g/cm³ 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³ 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³ 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³ 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³ buried at 3 m depth, with a width of 4 m, depth extent of 4 m, and a length of 500 m.

[0099] FIG. 9 is a plan view of the digital elevation model as measured from a LiDAR survey over the synthetic survey area.

[0100] 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.

[0101] 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³.

[0102] 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³.

[0103] 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.

[0104] 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.

[0105] 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.

[0106] 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.

[0107] 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.

[0108] 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.

[0109] 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.

[0110] 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.