Unmanned Airborne Magnetic and VLF Investigations : Effective Geophysical Methodology for the Near Future

Airborne geophysical investigations are now recognized as a powerful tool for geological-geophysical mapping, mineral prospecting, environmental assessments, ecological monitoring, etc. Currently, however, there are two main drawbacks to effective application of these investigations: 1) the difficulty of conducting geophysical surveys at low altitudes, 2) heightened danger for the aircraft crew, especially in regions with a rugged topography. Unmanned or socalled Remote Operated Vehicles (ROV) surveys are not bound by these limitations. The new unmanned generation of small and maneuvering vehicles can fly at levels of a few (even one) meters above the Earth’s surface, and thus follow the relief, while simultaneously making geophysical measurements. In addition, ROV geophysical investigations have extremely low exploitation costs. Finally, measurements of geophysical fields at different observation levels can provide new, unique geological-geophysical information. This chapter discusses future geophysical integration into ROV of measurements of magnetic and VLF electromagnetic fields. The use of GPS with improved wide-band Kalman filtering will be able to provide exact geodetic coordinates. A novel interpreting system for complex environments is presented that includes non-conventional methods for localizing targets in noisy backgrounds, filtering temporary variations from magnetic and VLF fields, eliminating terrain relief influence, quantitative analysis of the observed anomalies and their integrated examination. This system can be successfully applied at various scales for analysis of geophysical data obtained by ROVs to search for useful minerals, geological mapping, the resolution of many environmental problems, and geophysical monitoring of dangerous geological phenomena.


Introduction
The first airborne geophysical survey was carried out by Logachev in the former Soviet Union in 1936 [1].He used an inductor magnetometer to measure the magnetic field from an aircraft.This first air-magnetometer had a sensitivity of about 100 nT.By contrast, current aeromagnetic surveys in conventional vehicles can provide an accuracy of 5 nT and higher [2].Salem et al. [3] described an equivalent source approach that removes small environmental noise from the airborne measurements.Tezkan et al. [4] reported a recent pilot project that tested the feasibility of unmanned aircraft systems to carry out aeromagnetic surveys.Advances in reduced size and maneuverability of Remote Operated Vehicles (ROV) have led to an accuracy of 0.5 -1 nT (and higher) in surveys at low altitudes.The main advantages of ROVs are the absence of an aircraft crew (GPS controlled unmanned surveys can be carried out even at very low altitudes and in risky conditions) and comparatively low expenditures for geophysical observation projects.In addition, most methods of geophysical field analysis require knowledge of the field distribution at different levels over the Earth's surface.The ROV facilitates the collection of these data without having to use transformation methods.Finally, unmanned air geophysical surveys are extremely cheap (many tens of times less than conventional aircraft surveys).
A magnetic field is natural, and a Very Low Frequency (VLF) field can be considered to be quasi-natural [5] (given that perhaps several dozen VLF transmitters in the world radiate electromagnetic waves almost continuously).Combined aircraft magnetic/VLF measurements are potentially highly effective (e.g., [6][7][8][9][10][11]), but the highflying altitudes of conventional aircraft have limited the possibilities of using VLF.The ROV application can substantially reduce the distance between the airborne vehicle and the target while preserving the main advantages of conventional surveys.Thus the two main advantages of ROVs are their low cost and the fact that they are unmanned and hence do not put crew at risk.
ROV geophysical surveys at low altitudes can make important contributions to the search for minerals, geological mapping, various environmental (including archeological-geophysical) investigations, and military monitoring [4,[12][13][14][15][16][17].To increase the amount of data retrieved on such missions, new algorithms for the processing and interpretation of magnetic/VLF fields have been developed, especially for complex physical-geological conditions (mainly presented in [18]) including inclined magnetization (polarization), rugged terrain relief, unknown level of the normal field and a variety of anomalous sources.
Despite the progress in geophysical inertial navigation systems [19][20][21], etc., small ROVs flying at low altitudes are not sufficiently precise.The improved Kalman filtering with wide-band noise model presented here is designed to eliminate some of the obstacles to precise navigation.

Improved Kalman Filtering for Precise Vehicle Navigation
Kalman filtering (KF) is one of the most widely used methods for optimal estimation of random processes.The KF algorithm is the best linear real-time recursive filtering algorithm and is used to solve optimal filtration problems for non-stationary and multidimensional systems [22][23][24].However, the conventional Kalman filter is limited by the fact that KF uses white noise and that the useful signal and noise are independent.However, these limitations are a mathematical idealization that gives an approximate description of real noise.In fact, the actual noise is marked by a property that ensures the correlation of its values within a small time interval; i.e. if we denote this noise by , then ( ), 0 cov[ ( ), ( )] 0, where 0   is a small value.
The first attempts to introduce wide-band noise into the KF model to solve navigation problems in geophysics were carried out by Eppelbaum et al. [25].
The basic idea behind the method consists of modeling physical wide-band noise in the form of a distributed delay of the mathematical noise, i.e.

( ) ( ) ( )d t wt
where w is the white noise,  = const > 0, and  is the deterministic function.
with 0  s <  , i.e.  is the stationary wide-band noise.
In geophysical practice, forms of wideband noise can be expressed by autocovariance functions.Therefore, the problem involves modeling wideband noise in the form of (2), when the autocovariance function is given.
Let the wideband noise  have the autocovariance function (1), i.e., the function is known.Then from Equation (3) it follows that for   to be modeled in the form of Equation (2)  must satisfy the equation [26] (4)  Equation ( 4) is a convolution equation.It can be solved by direct and inverse Fourier or Laplace transforms.Application of the abovementioned approach can thus enhance the precision of GPS navigation in ROVs.

Elimination of Temporary Magnetic Variations
The accuracy of a high-precision ROV magnetic survey should not exceed the value of 1 nT.At the same time, the presence of both natural (basalts, diabases, gabbro, etc.) and artificial (iron and iron-containing) objects that have high magnetization may cause secondary variation effects which can be seen above this value.Conventional procedures for eliminating temporary (primary) magnetic variations are based on the trivial additive expression (e.g., [27]): where t is the time, x, y are the spatial coordinates, T (x, y, t) is the magnetic field recorded along a profile, 1 is the sum of "useful" anomalies and is the sum of the noise component caused by temporary variations.This formula only allows for simple subtracttion of the noise component from the observed magnetic field and cannot be used to calculate the effect of seconddary variations.
where T is the magnetic field, and b and c are the factors of a linear equation computed using the least-square method (LSM).
Since the secondary variation effect depends linearly on the intensity of the primary variations [28], the typical model of magnetic observations may be described in the following way [29]: , , , , , , It was shown that a b factor could be utilized to estimate the magnetization of the upper part of a geological section [18]: where are the field variations in time,  is the sum of the effects from the anomalous objects and geological inhomogeneities of the medium with regard to their dependence on the field variations.
where I is the magnetization of the upper part of geological section,  is the acute angle between the slope face and horizon, R is the slope length across the strike.
Given that  and R are known from the field observations and coefficient b is calculated using the LSM, the I value is easy to determine.Consequently, magnetic lowaltitude measurements along the inclined slopes may be used for quick estimation of the magnetization of the target medium.Eppelbaum [31] concluded that integrated ROV horizontal and inclined observations can serve to obtain the parameters of the medium magnetization even over a flat relief.
Measurements at the points of the profile made at different times t (2) make it possible to obtain a solvable set of algebraic equations so that the signal T can be extracted with sufficient accuracy.

Calculation of Terrain Relief Influence and Estimation of Magnetization of the Medium
Another disturbance factor is the influence of inclined relief on the magnetic measurements (the illustration here is limited to the case where the ROV follows the relief forms).Interestingly, this disturbance can help estimate the average magnetization of the medium.
The corrected field T corr may be formulated as: where T obser is the observed field.It is well known that in the case of direct magnetization, the field maxima correspond to ridges of the "magnetic" relief, while the minima correspond to its valleys [27,30].At the same time, the terrain relief effect is generally two-fold and can be attributed to the effect of the form and physical properties of the topographic bodies making up the relief and secondly to the effect of the slope of the observation line (profile), which causes variations in the distance from the field recording point to the target of excitation [18].The correlation technique below (see Equations ( 7)-( 9)) eliminates the effect of the first component.The second effect can be eliminated by using special expressions during the quantitative interpretation of magnetic anomalies [18].
Elimination of the topographic effect by the correlation technique allows practically complete smoothing out of the anomalies caused by the influence of uneven relief.
Along with the linear approximation of the relationship between the field and the height, approximations in the form of the square trinomial (parabolic equation) are appropriate here.This method can be used to improve the choice of the level for reducing the field to one plane.When computing a correlation field, the areas with the largest dispersion (compared to the averaging line), correspond to the profile intervals under which the disturbing bodies are situated.The presence of dispersion in itself is indicative of a hidden inhomogeneity in the section.An analytical approach [18] shows the possibility of obtaining the linear relation between the magnetic field and the height of observation to a typical element of rugged relief-a slope (an inclined ledge, or step).All the main relief forms can be approximated with the use of one or another combination of slopes (for each of these slopes, the correlation field can be constructed).Hence, as rough as it may be, there is a simple and effective method for eliminating the effect of magnetized rock relief.To apply this method, a correlation field is constructed between T (total magnetic field) and h (height of observation) values.For typical cases the correlation field is used to determine the terrain correction:

Initial Analysis of Magnetic Data
Visual analysis of magnetic anomalies is one of the most important steps in an investigation and all further stages may depend on the results of this initial decision-making.The standard principles of magnetic map analysis are presented in Table 1.However, for oblique magnetization (polarization), principles 3 -7 are not applicable.The shape of the anomalies is complicated by the horizontal component effect.This is why in the northern hemisphere, anomalies become asymmetric when the body's magnetization is parallel to the Earth's field.Thus, anomaly maxima of sublatitudinal-oriented bodies are appr Table 1.Common principles of magnetic map analysis (developed for vertical positive redundant magnetization) (after [18], with modifications).

No PRINCIPLES CHARACTERISTIC FEATURES OF ANOMALIES 1 Detection
Positive anomaly testifies to the presence of an anomalous object with increased magnetization (in the northern hemisphere).

Coulomb
Magnetic anomaly intensity directly depends on the excess magnetization of the object, its size and its depth.

Correspondence
The position in plane and shape (strike) of the anomaly coincides with the position in plane and shape (strike) of the anomalous object.
4 Maximum (consequence of the correspondence principle) Anomaly maximum is located close to the projection of the body's center (middle of the upper edge of the anomalous body) to the Earth's surface. 5 Gradient (consequence of the correspondence principle) The largest magnitude field gradients are observed close to projection of lateral boundaries of the body on the plane.

Symmetry
Plots of the field along directions across the isoline strike are symmetrical about the vertical straight line which traverses the maximum of the anomaly.

Limitation
The presence of minima on both sides of the higher maximum testifies to a dip-limited body.
shifted southward, whereas the minimum is located in the northern periphery of the anomaly.If the magnetization inclination with respect to the horizon is not large (in southern latitudes), the maximum is shifted more to the south, and the intensity of the minimum increases.Inflection points (maximum gradients) in the anomaly plots are displaced southward from the projections of lateral sides of the anomalous body on the plane.In the case of a depth-limited body, one of the minima (the southern one) may disappear along the maximum periphery.In the case of three-dimensional bodies of approximately isometric section, magnetic anomaly distortions occur not only due the inclination of the magnetization vector to the horizon plane, but also due to the difference in orientation of the horizontal magnetization projection with respect to the body's axes [32].Therefore, the analysis of graphs is not sufficient, and generally maps need to be analyzed as well.

Advanced Quantitative Interpretation of Magnetic Anomalies
One of the main interpretation problems is the quantitative determination of magnetic anomalies (inverse problem solution).The interpretation problem consists of a detailed description of the anomalous source using the measured magnetic field.However, the major principles of quantitative interpretation that were formalized for vertical magnetization do not work under conditions of oblique magnetization in low and central latitudes.Be-sides the influence of oblique magnetization, additional T anomaly distortions occur due to rugged (inclined) relief and unknown level of the normal magnetic field.This stage involves the application of rapid methods for quantitative interpretation of magnetic anomalies to construct an initial model of the geological section.Unlike certain conventional techniques [27,[33][34][35][36][37], the methods presented below are applicable in conditions of rugged terrain topography and arbitrary magnetization of objects where the normal field level is unknown [18].
In conditions of oblique magnetization, the "reduction to pole" calculation of "pseudogravimetric" anomalies [38] is often used.However, this procedure is only suitable when all the anomalous bodies in the area under study are magnetized parallel to the geomagnetic field, and simultaneously when the bodies have subvertical dipping.The magnetic fields can only be recalculated correctly when these restrictions are adhered to; this yields graphs that are symmetrical, and further interpretation using conventional methods can be carried out.Similar appro-aches based on the transformation of the observed magnetic field (for instance, analytic signal [39]), have the same limitations.
The methodology developed for complex environments [18] employs further modifications of the characteristic point method, the tangent method and the areal method, and utilizes the most commonly applied geometric models such as thin inclined bed (TIB), horizontal circular cylinder (HCC) and thick bed (TB) (a comprihensive explanation of these methods is given in [18]).These three geometric models, with different modifications, may be used for approximation and the corresponding quantitative interpretation of anomalies generated by various geological and environmental objects.
The procedures for interpreting anomalies generating the abovementioned models are presented in Tables 2-4.
The modern interpreting system developed for magnetic field analysis [18,40,41] includes the following components (besides the conventional ones): 1) Elimination of the secondary effect of time variations, 2) Calculation of the terrain relief influence and estimation of magnetization of the medium, 3) Application of an information-heuristic approach to geophysical field quailtative interpretation, 4) Inverse problem solution for complex environments (inclined relief, oblique magnetization and unknown level of the normal magnetic field), 5) 3-D modeling of magnetic (and gravity, if necessary) field.In the last 6) stage all kinds of information obtained in the previous stages (1)(2)(3)(4)(5) and conventional steps are integrated, and a final physical-geological model is developed.
The first stage of a ROV magnetic survey must include a priori investigation (generalization) of all available Table 2. Formulae for quantitative interpretation of magnetic anomalies over anomalous bodies approximated by thin seam and a horizontal circular cylinder using the improved characteristic point method (after [18], with modifications).

Parameters used for anomalies
Resulting from models Formulae to calculate parameters in terms of the anomalies resulting from models Parameters to be determined Thin bed Cylinder Thin bed Cylinder , , where cos 60 / 3 where 2 3 cos 2 sin cos / , where , where 3 3 2 cos 30 3 Indices "0" and "c" designate the thin bed and horizontal circular cylinder (HCC) models, respectively.Values h 0 and h c are the depths to upper edge of thin seam and center of the HCC, respectively.Parameter h designates measurements of magnetic field at different levels over the Earth's surface.magnetic properties of the targets and the surrounding medium.Figure 1 depicts a typical flow chart for aeromagnetic field analysis.
The following parameters are taken from the anomaly plot in the characteristic point method (Figure 2 In the tangent method four tangents are employed: two horizontal lines with respect to the anomaly extrema and two inclined lines passing through the points of the bend on the left-and right-hand branches of the anomaly plot.The following terms are taken from the plot (Figure 2): Table 3

. Formulae for quantitative interpretation of magnetic anomalies over anomalous bodies approximated by thin seam and horizontal circular cylinder using the improved tangent method (after [18], with modifications).
Parameters used for anomalies resulting from models Formulae to calculate parameters in terms of the anomalies resulting from models Parameters to be determined Thin bed Cylinder Thin bed Cylinder Generalized angle  where Parameters of magnetic moment M e , epicenter location and normal background level T backgr are determined from the values of θ and h obtained as in the characteristic points method by the formulae in Table 2.
Table 4. Formulae for quantitative interpretation of magnetic anomalies over anomalous bodies approximated by a thin bed and horizontal circular cylinder using the improved characteristic areal method [42].

Anomalous bodies Analytical expressions for calculation of characteristic areas Formulae for calculating M e and h
Thin seam where For a rough estimation of values h and M e the average values (without calculation of parameter Q 2 ) q h = 1.8, q M = 1.4 may be used.Then . Parameters of the magnetic moment M e , epicenter location and normal background level T backgr are determined from the values of θ and h as in the characteristic points method by the formulae in Table 2.
of an inclined tangent with horizontal tangents on one branch; d 4 = the same on the other branch (d 3 is selected from the plot branch with conjugated extrema), d 6 = interval between d 3 and d 4 , d 8 = difference in abscissae between point of intersection of left inclined tangent with lower horizontal tangent and inflection point on left branch (d 6 and d 8 are used for the model of thick bed only).A detailed description of magnetic field examina-tion of a thick bed model is given in [18].
The areal method is based on calculation of separate areas limited by the anomalous curve, horizontal line and two vertical lines crossing some singular points at the anomalous curve.
When anomalies are observed on an inclined profile, the obtained parameters characterize a fictitious body.
he transition from fictitious body parameters to those of T Copyright © 2011 SciRes.POS the real body is done by applying the following expressions (the subscript "r" stands for a parameter of the real body): where h is the depth of the body upper edge occurrence (or HCC center), x o is the shifting of the anomaly maximum from the projection of the center of the disturbing body to the Earth's surface (caused by oblique magnetization), and  o is the angle of the terrain relief inclination ( o > 0 when the inclination is toward the positive direc-tion of the x-axis).Figure 3 illustrates the application of these methodsimproved tangents and characteristic points (Figure 3(a)) and characteristic areal method (Figure 3(b))-for the airborne magnetic field analysis in the Big Somalit (southern slope of the Greater Caucasus).

A Scheme for VLF Data Examination
The most effective scheme for VLF data observations is when the magnetic and electric components of the electromagnetic VLF field are measured.The VLF equipment developed in the former USSR and Russia e.g., [43,44] can be used to obtain the values of magnetic compo- nents on coordinate axes: H x , H y and H z , as well as changes of radio-bearing  to the VLF transmitter (as well as other parameters).The use of magnetic VLF components in many cases is more efficient because of the comparatively low level of the various kinds of noise [45].Observation parameters such as tilt angle and ellipticity are not discussed in this article.
An important peculiarity of a VLF ROV survey at low altitudes is the slow decrease in the VLF signal with increasing observation height over the Earth's surface [43,46].

Elimination of Temporary VLF Variations
The problem of eliminating time variations in the VLF method is crucially important since the intensity of these variations is often compatible with the intensity of the averaged VLF signal during the time of the field survey.The following approaches have been used to overcome this problem.
In early VLF methodology, an approximate technique was used based on field intensity measurements at a control point (CP) before and after fieldwork [47].The corresponding corrections were introduced by interpolation assuming that the field intensity change was linear.However, practical application of this method led to con- siderable error.
A modification of this method is described in [48].The difference is that the time period of supposedly linear variations is assumed to be 1 hour.The area is rejected if the level of variation for 1 hour exceeds 20%.This method, however, suffers from following disadvantages.Experience in field explorations in various regions indicates that the intensity changes of the VLF fields under study (even over the course of one hour) often cannot be approximated by a straight line with sufficient accuracy.In addition, the amount of variation during this period may exceed 20%.It is also worth noting that this method neglects the effect of variation noise intensities on useful anomalies.It is obvious, however, that time variations in VLF field intensity affect the radiowave energy passing from the air into the ground.The intensity of the secondary field tends to increase or decrease, respectively (especially in the presence of ano-malous objects of increased or decreased conductivity), and the variations are non-uniform.When disregarded, this fact can distort the interpretation of the results.
The Scintrex Company developed an automatic attachment to a VLF receiver.Variations are eliminated by obtaining a synchronous ratio of the vertical component H z of the VLF magnetic field to the total horizontal component H  .However, it was found by experimental field exploration that the H z /H  ratio may vary even for observations at the same point.Quantitative interpretation of the H z /H  curve presents certain difficulties.Another approach described by Scintrex in a recent guide for using the VLF equipment, consists of the selection of the frequency and time interval that is the most stable temporally for the target region.Vallee et al. [49] suggest a similar method.However, this approach would benefit from improvement.
The wavefront is thus assumed to be a plane whereas the primary field intensity is taken to be constant within the area of the detailed geophysical investigation (25 -100 km 2 ).The large distance from the utilized VLF transmitters (2,000 to 10,000 km) accounts for the uniformity of the primary field [46,50] both at CP and in the investigated area.
The following additive model of a geophysical field is applicable in practice for investigation by the VLF method:  , where F j is the field observed along the profile,  is the sum of effects from the anomaly-forming objects and geological inhomogeneities of the section (taking the dependence on field variation into account), represents the noise of field variations in time (time dependence is omitted in Equation (11) for simplicity).

 ,
It is known that for the VLF field the conduction currents dominate the displacement currents [43,51], which is the major condition for quasi-stationarity [52].The proportionality of the secondary electromagnetic fields (both magnetic and electric) to the primary fields has been confirmed in a number of publications e.g.[53,54].Hence given the physical effect of the subvertical radiowave propagation in the ground [43,55] and taking into account Equation (11), the following simplified model of VLF observations can be used in practice [56]: where H o (t) is the observation at the control point; H j (t) is the observation on the profile, H(t) is the primary field intensity, B and C are coefficients reflecting electromagnetic properties of the medium; indices "j" and "o" mark the observation point on the profile and CP, respectively.Considering that B o = const, this parameter, which is a basic value for the given area, can be assumed for convenience to be equal to zero.Solving Equations (12) in the parametric form, we obtain the values H clear cleared from the variations for each profile point: clear ( ) where H is the certain averaged value of the field.Thus, the proposed scheme of eliminating variations in the VLF method is based on synchronous recording of the observations along the profile and variations at CP. Automatic variation recordings present no difficulties and can be carried out with a discretization interval of 3 -5 seconds.This method substantially eliminates observation distortions caused by field variations with time (both over the course of a day and at different days in a survey) by reducing observation results to a common level [57].
It should be noted that Bozzo et al. [58] later noted the necessity of normalization of field VLF measurements to the primary field observations.

VLF: Calculation of Rugged Relief Influence
Attempts have been made by many investigators to reduce the terrain relief effect.An analytical approach was suggested by Tarkhov [47]; however, the analytical calculations proved to be too cumbersome even for a model of a uniform slope made up of homogeneous rocks and, above all, were at variance with the experimental results.Sedelnikov [59] analyzed in detail the theoretical problem of electric field distortions caused by a single mountain and a mountain range but suggested no method for eliminating topographic effects.
Karous's [60] publication about topographic influences on the vertical component of the VLF magnetic field (H z ) is clearly of interest, but no numerical results were reported.Eberle [61] proposed simplified formulas to calculate relief corrections for the observed electromagnetic fields; their application is however limited by a number of conditions, e.g., the requirement that the disturbing object and the observation profile should be in the most favorable position for excitation by the primary VLF magnetic field.This requirement is seldom obtained or difficult to meet.
The second approach is based on a physical (analog) simulation of electromagnetic fields.Gordeyev et al. [43] constructed models which describe quantitatively the distribution of the VLF electromagnetic field over different terrain relief forms.Baker and Myers [62] suggested an electromagnetic field "reduction" method based on their model investigations.This method, however, has  the same 1imitations as found in [61].
Gordeyev [63] proposed a third approach.It is based on using a graticule (chart) method to reduce the topographic effect on the total horizontal magnetic field component of remote transmitters.It is difficult to apply at greater relief dips and can, in some cases, impair the interpretation.This limitation is traceable to the "build-up effect" of a rectangular graticule in a complicated (relative to the observation profile) topography, which in addition is often characterized by different electric properties.
It should be noted that the effect of the Earth's surface topography on radiowave propagation is a fairly important problem in radiophysics and has been investigated by quite a number of researchers.An overview of the literature suggests that this problem has yet to be adequately resolved.A numerical calculation of radiowave scattering is a complicated electrodynamic problem even for the simple case of diffraction by a homogeneous wedge displaying finite conductivity [64,65].The difficulties increase progressively for statistically inhomogeneous (i.e.real) surfaces.Bass and Fuks [66] emphasized that under natural conditions there is no need to find the detailed structure of a scattered field and that it is sufficient to ascertain some reflected-signal parameters averaged for the whole class of surfaces and targets.In general, researchers have apparently come to the conclusion that wave scattering is determined largely by coordinate functions and represents a combination of wave diffraction on an arbitrary surface and an application of the theory of probability and mathematical statistics.
A review of the literature and results of VLF surveys lead to the following general rules: 1) The plots of H x (the x-axis is aligned with the observation profile) and the total horizontal magnetic field component H  substantially copy the surface topography.There is a direct correlation between the increase or decrease in the recorded field intensities and positive or negative landforms, respectively.2) The intensity of the recorded VLF magnetic field is directly proportional to the relative elevation of observation points.The topographic anomalies caused by local relief forms are similar to those from conductive bodies but have a smoother appearance.3) When the angle between the radiowave arrival direction and that of a relief element is close to 0˚, the relief effect is relatively large and when this angle is close to 90˚, it is relatively small, but in any case the relief effect is always present.Since all observations are made in the "distant" ("wave") zone, the bearing changes negligibly along the profile and the topographic anomaly is quasi-independent of the bearing change.
4) VLF field intensity is mainly dependent on the resistivity of the rocks in the upper portion of a geological section.
All the abovementioned points make it possible to apply the correlation method developed in magnetic prospecting (see section "Calculation of terrain relief influence and estimation of magnetization of medium") to the VLF method [45].An example of the correlation approach is presented in Figure 4.This example illustrates the application of the correlation technique on a portion of the Katekh pyrite-polymetallic deposit in the southern slope of the Greater Caucasus.A frequency of 16.0 kHz was used in the investigation (the transmitter was in Rugby, Great Britain), and the profile azimuth was 60˚.The ribbon-like band of the pyrite-polymetallic ore was located at a depth of about 60 -80 m in the sandy argillaceous series of the Upper Aalenian.The average resistivity of the nearly uniform surrounding medium was in the range of 700 to 900 Ohmm; the resistivity of the massive ore body was a fraction of 1 Ohmm.The skin depth in the surrounding medium was about 110 m.Since a topographic anomaly was superimposed on a small signal from a relatively deep-seated anomalous object, it was difficult to detect the latter.
Two different correlations were made for the southwestern and the northeastern slope, yielding correlation coefficients (between the relief heights and VLF intensity) of r 1 = 0.988 and r 2 = 0.
. After removing the relief effect, a positive anomaly can clearly be seen on the northeastern slope in the corr plot.This anomaly may be due to the edge effect of a deep ore body, whereas the minimum on the southwestern slope may correspond to the opposite edge of the subhorizontal tabular ore body, since the anomaly is not very large in this case [45].

VLF: A Brief Review of Available Methods of Quantitative Analysis
An overview of the literature indicates that there are practically no reliable or rapid techniques of quantitative interpretation in the VLF method.The methods suggested by Tarkhov [47] are semiquantitative.The procedure developed by Fraser [67] is useful primarily for revealing anomalous objects.
The techniques presented by Gordeyev and Sedelnikov [68], based on calculating the anomaly extrema ratio, require knowing the zero line (or a normal field).The wrong choice of the zero line entails substantial errors in determining the anomalous object's quantitative parameters.Moreover, these techniques are not intended for interpretation in conditions of rugged terrain relief.Dmitriyev et al. [53] solved the direct problem for a number of electromagnetic methods (including the VLF method).They provide a numerical analysis of the anomalies due to bed-like bodies as a function of their electric properties, dimensions and position with respect to the Earth's surface and the observation profile.However, computational difficulties restrict the range of models to simple ones.
Zhdanov and Keller [69] noted that the finite differences method is the most effective method of electro-magnetic field mathematical simulation.However, the expediency of its application depends on the size of the design area and the desired accuracy of the computations.For a large area of calculations or when the simulation problem needs higher accuracy, the computations exceed all possible limits even for high-performance modern computers.
Khesin et al. [18] pointed out that some approaches developed for magnetotelluric sounding [69,70] could be employed in the VLF method.Basokur and Candansayar [71] suggested using the same methodology (however, these authors concluded that this interpretation could only be qualitative).However, the application of these procedures is usually limited by the extreme variations in electric properties in the (upper) part of geological section.
Gordeyev and Sedelnikov [68] and Gordeyev et al. [43] characterized the physical (analog) simulation of VLF curves by models of a vertical inclined thin bed.The accuracy of the physical simulation and that of the numerical computations (for simple models) is at present approximately the same, and is roughly 8% -10%.
Numerous other works, in many respects similar to those mentioned above, have dealt with the anomaly interpretation in the VLF method.For instance, Karous and Hjelt [72] and Olsson [73] reported the calculation of plots for vertical and horizontal components of the magnetic field of VLF using a model of an inclined thin bed with different dip angles.The depth of occurrence of the upper edge of the bed varies as well.Olsson [73], Tesmull and Crossley [74] and Poddar [75] computed components of the VLF fields for different types of geoelectric sections.Miecznik [76] calculated the magnetic components of the VLF field caused by another model of an anomalous object-that of a conductive cylinder placed in a homogeneous medium (the cases of both E-polari- zation and H-polarization are discussed).
The characteristic VLF diagrams developed by Sinha [77] can provide useful information about vertical and inclined conductive beds in simple geological-geophysical conditions.The VLF first derivative method proposed by Djeddi et al. [78] is based on perspective (the authors suggest eliminating the topographic effect of peculiarities in the interpretation process).
Olsson [55] solved the direct problem for an ideal conductive half-plane with different dip angles, when the half-plane is covered by variable-depth loose deposits.In addition, he proposed techniques for a simplified interpretation of VLF data based on the utilization of the extremum points in the anomaly plot.The synthetic models presented Pedersen and Oskooi [79] provide a better understanding of the resolving power of the VLF data.The same authors propose employing a tensor variant of VLF measurements.Monteiro et al. [80] described a quantitative interpretation of VLF-EM data using 2-D models (the authors developed a 2-D regularized modeling of VLF-EM data based on a forward solution using the finite-element method) that yields practical results about the Earth's average resistivity in the survey area.Kaikkonen and Sharma [81] performed a 2-D simulation of noise-free and noisy media with single-body and complex models.The authors highlight the need to use reliable a priori information before running such a simulation.

VLF Field: Possible Application of Procedures Developed for Potential Field Analysis
The main equations used in the theory of alternating electromagnetic field are Maxwell's Equations: Here H and B are the magnetic field vectors (H stands for intensity and B for induction); E and D are the electric field vectors (E stands for intensity and D for induction); q is the charge density; j is the conduction current density.
Equation (14a) can be written as [52]: where C is the total current density; j cond and j disp are the conduction and displacement current densities, respectively.It has been found that at frequencies of 10 to 30 kHz, j disp  0 [51].Therefore, rot H = j cond (16) or where is the electric conductivity.
An EM field can be considered quasi-stationary if it satisfies the three conditions of quasi-stationarity.Landau and Lifshitz [82] described the physical meaning of these conditions.Alpin et al. [52] formulated these conditions as applied to geophysics.These requirements are as follows: 1) Slow field variation.
3) To avoid an appreciable lag in the magnetic field variation, the region including magnetic field generating currents and observation points must not be too large.Let us consider these conditions in more detail.
1) The study of time variations of the VLF fields has shown that sharp changes in intensity are not characteristic of these fields.Nonetheless, variations that use a special procedure in fieldwork should be eliminated (see section "Elimination of temporary VLF variations").
2) This condition follows from Equation ( 16), since div rot H = 0 and, therefore, div j = 0. Almost closed currents do not violate this condition.
3) This condition can be formulated in the following way: the length of an electromagnetic wave in the ground should essentially exceed the length of the investigated objects.This condition, as a rule, is fulfilled in the prospecting of mineral deposits, geological mapping and object analysis in engineering and military geophysics.
It is known that the fundamental solution of the Laplace Equation in the 2-D case is the following function: where R MP is the distance between points M and P (M is the observation point and P is the point of the body).
The fundamental solution of the Helmholtz wave equation is a Hankel function of the first kind with zero order (1) 0 H :

  
(1) 0 0 , 4 where k is the wave number and i is the imaginary unit.The function (19) is the analog of the function (18) for a Laplace equation.
Using the above correspondence, Hänl et al. [83] and Dmitriyev [54] proved the relationship between Green's solutions of the Laplace and Helmholtz equations.This allowed them to show how the results obtained in potential theory could be extended to the Helmholtz equation.
Dmitriyev et al. [53] indicated that the association of the singular points resulting from the VLF field anomaly plots with geometrical parameters of the bed is analogous to the well-known behavior of singular points in potential field anomalies.Zhdanov [84] carried out a theoretical investigation of the application of a set of Cauchy-type integral analogs to the problems of electromagnetism.It was suggested that the methods developed in potential theory for analytical continuation, separation and quantitative interpretation could be applied to quasi-stationary electromagnetic anomalies.
The quasi-stationarity of the electromagnetic field follows from condition 1).

POS
An overview of publications [43,75,[85][86][87] indicates that a conductive thin inclined bed (TIB) is the most common model for the VLF technique.The analytical expressions of the fields for this model (the 2D case) are presented in Table 5.
In this table Z v and X v are the vertical and horizontal magnetic field components for vertical magnetization, respectively; H x and H z are the horizontal and vertical VLF magnetic field components, respectively; H o is the VLF primary field intensity; k is the coefficient reflecting geometry and conductivity of the bed.
Evidently, Z v is proportional to H x , whereas X v is proportional to H z .Similar results were obtained for a HCC model [44].
Thus, the plots of magnetic fields of VLF transmitters can be interpreted by special methods elaborated in potential field theory (in particular, in magnetic prospecting, for conditions of rugged terrain relief, oblique magnetization (polarization) and an unknown level of the normal field).The E-polarization vector, in the first approximation, is the analog of the magnetization vector [44].
In this context, the special techniques presented in Tables 2-4 are of practical interest.In the VLF method can be applied the interpretation techniques developed in magnetic prospecting for inclined thin bed and horizontal circular cylinder models.It should be noted that the total horizontal component of the VLF magnetic field The essential distinctions reflecting the quasi-potential nature of the VLF electromagnetic field are as follows.In magnetic prospecting the induced magnetization vector of an inclined bed is approximately parallel to the geomagnetic field vector, irrespective of the bed dip direction, where the magnetic susceptibility does not exceed 0.1 SI unit.In electric prospecting by the VLF method under E-polarization for highly conductive objects, the equivalent vector of polarization that causes the anomalies of H x and H z , approaches the body axis, with a slight deviation toward the vertical for gently sloping bodies [43,88].This makes it possible to utilize the generalized angle θ derived from these methods and to represent the difference in the inclination angles for the bed and the polarization vector to estimate the bed dip angle by the following empirical formulas [44]: 1) for the H x anomaly 3 90 2) for the H z anomaly 3 180 where  is the bed inclination.
For the case of observations on a sloping relief, the angle  is calculated as follows: 1) for the H x anomaly   2) for the H z anomaly   ( o is derived in Equation ( 10)).
It should be noted that the position of the upper edge could be slightly shifted from the real upper edge downward along the bed dip.This is because the linear currents are concentrated in the upper portion of the conductive object.This portion may be situated below this object's upper edge.
An actual example of VLF data analysis (land survey) is presented in Figure 5.It should be noted that the anomalous magnetic X-component computed for the same model is very similar to the observed VLF H z curve (their similarity follows from Table 5).The figure clearly shows that the results fit well with the geological data.

Revealing Buried Ring (Circular) Structures
One of key problems of airborne magnetic and VLF prospecting is the delineation of typical buried objects of various size and origin.Ways of revealing such buried objects with an a-priori shape are illustrated below through examples of identification of ring (circular) buried targets.
Ring structures (RS) are omnipresent in the Earth's environment (Figure 6).RS are generally classified as Terrestrial and Extraterrestrial (natural RS) and Archaeological and Military/Engineering (artificial RS).It should be noted that the dimensions of natural RS are often compatible with the size of artificial RS (excluding well-studied RS measuring tens or hundreds of km in diameter) [42].Thus problems sometimes arise not only of RS identification, but also of correct classification.Greater accuracy can be achieved by collecting more comprehensive geological (geochemical), geophysical, archaeological and engineering (military) data.In the flow chart below (Figure 6) natural RS are classified using Khain's [89] scheme and the author's supplements and extensions [90].Archaeological RS are presented as buried caves, roads and circus rings, remains of fortresses and towers, tailings of detritus and religious edifices.Detection should be coordinated with conventional land geophysics.Military RS are displayed as rocket shafts, command posts and defensive installations (obviously, classes of these RS may be significantly extended).In the latter case the ROV application has clear cut advantages over other identification methods.

Methodology of RS Delineation
The first problem in RS analysis is to differentiate the targets from background noise.Very frequently it is difficult to single out RS in complex geological-geophysical environments, especially when the typical rectangular network of geophysical observations is used.In addition, anomalies caused by RS can be distinguished from geophysical field features by their location on the periphery of anomalous bodies.For this purpose a special method has been devised to distinguish concentric structures by computing the sums of the horizontal gradients of geophysical fields and their differences [18].It should be noted that ROV applications make it possible to observe geophysical fields at different levels and compute the difference in horizontal gradients between them, which might be employed as additional searching indicator.
The apparent graticule radii drawn at an interval of 45˚ determine the horizontal gradients (Figure 7(c)).The sum of the gradients should be higher in the presence of circular features and other signals should remain constant.
Here the correlation of the sum of the gradients (or the average gradient) for a circle with a radius R n and a ring external to this circle limited by R n and R n+1 radii makes it possible to determine whether the circular feature reflects a centric or ring structure (Figure 7(d)).The sum of gradients inside the circle tends to zero in the absence of a centric texture [18].An application of this method is depicted on a model of an inclined circular cylinder magnetized along its dip (Figures 7(a) and (b)).

Multimodel Approach to Geophysical Data Analysis
The multimodel approach to geophysical data analysis is described in the example below of detailed magnetic prospecting.A quantitative interpretation of magnetic anomalies (and many other geophysical methods) was traditionally oriented to a single model for buried object identification.When several hypotheses relating to the parameters of the body causing the disturbance (i.e., the buried object) were made, usually only one model was selected that roughly presented the object in the domain  x of k-dimensional space of the physical-environmental factors.As a rule however, human, industrial and geodynamic activity as well as various geological/environmental processes affect many environmental features.Additional sources of noise that affect interpretation include rugged terrain relief, oblique polarization of geological objects and a heterogeneous host medium.Thus the response function  i -the geophysical field-may ambiguously represent an archaeological target.To overcome this issue, domain  x can be divided into several subdomains  1 ,  2 , ... ,  m and in each a single model will dominate [92].In such a way m physical-environmental models of the same target can be developed, where each is corrected for separate subdomains  1 ,  2 , ... ,  m .
For instance, remains of an ancient wall in specific geological conditions can be defined as a thick bed (magnetic field), a sphere (gravity field) and a horizontal plate (VLF or resistivity).The multimodel approach can also be applied at varying levels of geophysical field recording.As a result, different explanatory models contribute to the process of quantitative interpretation.Integrating several response functions  I , leads to a more accurate and reliable physical-environmental model of the buried target.
The simple model presented in Figure 8 shows how wo different interpretations of the same ancient ruin can t  Combining these two models (we have two response functions  1 and  2 from the subdo-mains  1 and  2 ) yields a combined generalized model of the anomalous body [17].Similar models could be developed for various scales and different geophysical methods.
In general combining different geophysical methods with their multilevel (including different depths of electrodes or geophone grounding) observations [16,17] produces a large number of combinations, which may be used for extracting additional useful information from the ata.d

Integrated Analysis of Magnetic and VLF Fields
Integrated analysis of magnetic and VLF fields can be carried out using both informational and wavelet approaches

Informational Approach
Since the solution of geological problems calls for inte-gration of different geophysical methods, the abovementioned features (geophysical observations, petrophysicalmorphological characteristics and interpretation process) should be expressed in the same terms.Only then can specific concrete problems for each method be formulated.At a specified level of accuracy of geophysical meth- Copyright © 2011 SciRes.POS ods and object characteristics, the integrated interpretation technique is of primary importance for integrating geophysical methods.The same is true when defining requirements and preconditions for data collection and the assessment of the findings from each method.This information can be obtained from physical sources corresponding to various classes of geological (environmental or technical) objects.The integrated interpretation, which leads to class intersection, makes it possible to single out one or several geological targets, which form a solution to the geological objective.In some cases only the comparison of fields measured by various geophysical methods leads to an "informational leap" [93], which can then single out the target objects.In some cases, only an estimate of the amount of corresponding information contained in observations can substantiate the presence of an object in a given class.These estimates are especially valuable for integrated interpretation, as they provide a means to operate on data obtained through various methods which are now expressed in a common (informational) unit.The target is determined in terms of the maximum information obtained from the geophysical methods set.
The language of information theory [94] is the most suitable to express the essence of geological and geophysical investigations, i.e. the steps in the acquisition and processing of information.Therefore, in spite of the probabilistic calculus of the amount of information, the informational approach to this problem introduces deterministic features.Here, the advantage of the logical nature of a number of interpretation methods based on information-statistical (logical-statistical) approach becomes evident.Furthermore, this approach includes such probabilistic methods as correlation, filtering and statistical decision-making, which are widely accepted for data processing and interpretation.
Geophysical data processing is mainly intended to reduce and eliminate different kinds of noise.The main task of qualitative interpretation is to single out an object of the desired class, while that of the quantitative interpretation is to determine and refine parameters of the object.The solu-tions of these problems are closely interwoven and based on the model conceptions of the interpreter.
At each point the amount of information J i due to the application of the i-th geophysical method may be expressed as [93]: or where P j is the relative frequency of the j-th interval of the i-th indicator on the histogram of its distribution, U i and U i are the amplitude and the error of this indicator's determination, respectively.Even such a simple calculation of the amount of information may be effective one when operating with fields that differ physically.
After summing up the information elements, which can indicate a priori that the object of the desired class is , can be computed.When a large amount of information is contained in the data of only one or two methods, an integrated criterion J int can be computed (it is only one of several types of information integration).This choice depends on the number of significant indicators and whether their relative influence is reduced [93]: where J p is determined using pairwise combinations of the results of n methods used (I 1 and I 2 are pairwise geophysical observations from the total number n of applied geophysical methods): To avoid missing the detection of deeply buried objects, it is worthwhile in some cases to use the relative frequencies of average values or average field estimates on a sliding averaging interval instead of the P j and U i values, respectively.
Along with a simplified version based on summing up informants obtained from Equation ( 23), the method also exists as an "Integration" program.The latter can compute the sums of informants by Equation (22) and the corresponding parameter J integr (Equation ( 24)).Unlike other methods, this index is capable of revealing the objects characterized by the maximum number of indicators of different intensities and, at the same time, avoids missing an object which for some reason did not manifest itself by any indicator.The combination of indices serves to make interpretative conclusions.
In practice J i is usually replaced by the relative amount of information, which is also known as the coefficient of informativity The value of i J determines the information obtained when the result of U j falls into the x j interval of the histogram with an equal probability of falling into any of the R intervals.According to [95], this is equal to average (complete) information, contained in the results achieved by measuring with a single method, log .
Application of K i takes into account differences in the ranges of different fields.However, the application of expressions ( 22) and ( 24)-( 26) may be not effective when there are an insufficient number of geophysical observations.

Wavelet Approach
About 20 years ago a new branch of computational mathematics-Wavelet Analysis (WA)-was developed.WA has had a significant impact on the cardinal problems of signal processing such as denoising of signals, enhancement of signals and differentiating signals with closely related characteristics [96][97][98][99][100][101].
Researchers typically need to deal with the problem of extraction of essential features from available data contaminated by random noise and a non-relevant background.If the essential structure of a signal consists of several sine waves, it can be represented trigonometrically (Fourier analysis).In this case the signal can be compared with a set of sinusoids and the consistent ones can be extracted.An indicator of the presence of a wave in a signal   f t is the Fourier coefficient Wavelet analysis provides a rich library of available and rapid, computationally efficient procedures to represent signals and select relevant waveforms.The basic assumption behind wavelet analysis is that the essential structure of a signal consists of a small number of various waveforms.The best way to determine this structure is to represent the signal by a set of basic elements containing waveforms coherent with the signal.For structures of the signal coherent with the basis function, large coefficients are attributed to a few basic waveforms, whereas small coefficients are expected for the noise and structures unrelated to the basic waveforms.
Wavelets are a family of functions ranging from functions of arbitrary smoothness to fractals.
The wavelet procedure takes place in two steps.Step 1 is decomposition, i.e. breaking up the signal to obtain the wavelet coefficients, and the second is a reconstruction, which consists of re-assembling the signal from coeffi- is the continuous wavelet transform.

POS
The wavelet approach has been successfully tested on numerous models of karst terrains for which gravity and magnetic fields as well as GPR images were computed [101,102].We chose to characterize the geometric events in the merged data curves by a representation of the curves in the Fourier domain.The assumption is that each of the A (absence of object) and P (presence of object) datasets is characterized by a number of dominating frequency bands.A perfect tool to reveal these characteristic frequency bands is provided by wavelet packet (WP) analysis [99].Once implemented, the wavelet packet transform of a signal yields highly redundant partitions of the frequency domain.As a result, geophysicists can automatically detect the karst occupied areas and areas not containing this dangerous geological phenomenon on the basis of self-testing wavelet methodology.

Conclusions
Geophysical surveys using the latest generation of remote operated vehicles (ROV) have provided the impetus for mobile and detailed geophysical mapping of large territories including mountainous regions, swampy areas, areas of dense natural growth, military surveillance zones, and even archaeological sites.ROV lanes may be projected to levels previously inaccessible to geophysical aircraft to a meter above the Earth's surface.Utilization of wide-band noise in the improved Kalman filtering model can significantly increase the accuracy of ROV navigation through the use of GPS.The advanced magnetic data methodology under complex environments (oblique magnetization, rugged relief and unknown level of the normal magnetic field) presented here can be effectively applied not only to unmanned airborne magnetic survey data, but also, with some modifications, to ROV-derived VLF measurements.In addition, a new approach to removing time variations in the VLF method was suggested.A design for the delineation of buried ring structures (as one of the possible search targets), and the application of a multimodel approach (on the example of multilevel measurements) was discussed.An integrated analysis of magnetic and VLF data can be carried out by applying both informational and wavelet approaches.Thus multilevel integrated ROV magnetic -VLF surveys will doubtless become a powerful geophysical tool in the very near future.

Figure 1 .
Figure 1.A generalized flow chart for airborne magnetic field analysis.

Figure 2 .
Figure 2. Position of characteristic points and tangents at T anomaly due to an obliquely magnetized thin bed.
as the H x component, since the contribution of the H y component is usually relatively small.

Figure 5 .
Figure 5. Quantitative analysis of the vertical magnetic component H z of the electromagnetic VLF field in the copper-nickel deposit (Kola Peninsula, Russia) using the improved tangent and characteristic point methods (observed curve and geological section after [43]). 1) copper-nickel ore body, 2) location of drilled boreholes, 3) results of quantitative analysis: cross indicates position of the upper edge of the ore body and arrow shows its dip direction.
be obtained by conducting a magnetic survey at two different levels (0.1 and 3.0 m, respectively).The survey at the 0.1 m level generates a typical thick bed (TKB) model (Figure8(a)) and at the 3.0 m level, observations of the anomalous body may be interpreted as a horizontal circular cylinder HCC (Figure 8(b)).The results of the TKB model interpretation were used to determine the center of the upper edge of the anomalous body (Figure 8(a)) and the HCC model for localization of the center of HCC (Figure 8(b)).

Figure 7 .
Figure 7. Singling out a concentric structure (after [18]): (a) model field of an inclined circular cylinder calculated along a profile, (b) model field complicated by random noise (in plane), (c) apparent graticule for concentric structure selection, (d) singling out a model body by summing horizontal gradients of the field within apparent circular graticule zones.Isogams of a model field (b) and the sum of its gradients (d) in conventional units: 1) positive, 2) zero, 3) negative; cylinder edge projection: 4) upper, 5) lower; 6) contour of the portion shown in (b).

Figure 8 .
Figure 8. Realization of two-level observations with two different interpretation models: (a) model of a thick bed, (b) model of the horizontal circular cylinder (HCC).Magnetization is denoted as J.
are many variants of WA.In the so-called continuous WA the signal f(t) is tested for the presence of waveforms.Here, a is a scaling parameter (dilation), and b determines the location of the wavelet t