Improved techniques in data analysis and interpretation of potential fields: examples of application in volcanic and seismically active areas

Geopotential data may be interpreted by many different techniques, depending on the nature of the mathematical equations correlating specific unknown ground parameters to the measured data set. The investigation based on the study of the gravity and magnetic anomaly fields represents one of the most important geophysical approaches in the earth sciences. It has now evolved aimed both at improving of known methods and testing other new and reliable techniques. This paper outlines a general framework for several applications of recent techniques in the study of the potential methods for the earth sciences. Most of them are here described and significant case histories are shown to illustrate their reliability on active seismic and volcanic areas.


Introduction
Gravity and magnetism are classical branches of the earth physics.Potential field interpretation represents one of the oldest and most suitable approaches followed by geophysical researchers.This is also due to the high cost-benefit ratio which characterizes the application of these methods and by the relatively simple logistic management of a gravity and magnetic survey.
The amount of computation often required by numerical algorithms has long prevented the full use of many theoretical efficient methods of analysis and interpretation of the potential fields.This is the case of several methods involving computation of spatial derivatives.Many of these interpretation techniques fall in the class of the enhancement methods (Blakely, 1996).
In the following, improvements to well established enhancement methods (downward continuation, vertical and horizontal derivatives, Euler deconvolution and 3D boundary analysis) are described and illustrated by applications to Italian seismic and volcanic active areas.

Enhancement methods
The advantages of using the vertical derivative of gravity were first recognized by Evjen (1936).At the beginning of 80's, rather different directions were suggested to develop techniques useful to geologic interpretation of gravity and magnetic maps.The horizontal derivative method (Cordell and Grauch, 1985) is the best example of such techniques.By this technique, it is possible to locate the horizontal position of the density or magnetization boundaries.Afterwards, this method and others were implemented to give also an estimation of the depth to the sources.They assume particular source-models for which the depth estimates are valid, like the Werner deconvolution (Hartman et al., 1971), the analytic signal (Nabighian, 1984;Roest et al., 1992) and the enhanced analytic signal (Hsu et al., 1996(Hsu et al., , 1998)).
Information on the depth and other geometrical source parameters are provided by the method known as «Euler Deconvolution».It was first proposed by Thompson (1982) and extended to the analysis of maps by Reid et al. (1990).Thurston and Smith (1997), Smith et al. (1998) and Thurston et al. (1999) by their Source Parameter Imaging method obtained essentially the same information as the Euler deconvolution by computing the so-called 'local wavenumber'.
Finally, another rapidly developing family of methods is that based on the Continuous Wavelet Transform (CWT; see Moreau, 1995;Hornby et al., 1999, Sailhac et al., 2000).Using a particular class of wavelets, based on the Poisson Semigroup kernel, these authors explored the information contained in the continued and derived field.In such a way 2D magnetic problems may be solved for both the source depth and the homogeneity degree that is in turn related to the source shape.The determination of depth and other source parameters are here obtained from the knowledge of the field at several altitudes (or scales), similarly to what was already achieved by other authors (Paul et al., 1966;McGrath, 1991;Fedi and Rapolla, 1999).

Stable downward continuation by Integrated
Second Vertical Derivative (ISVD) method: an application to the gravity field of the Campanian Plain (Italy)

Computational approaches
The shape of a gravity and magnetic anomaly is intrinsically smoothed with respect to that of its source.A physical way to improve resolution of potential fields is to calculate the field to a level close to the sources by means of the downward continuation.
Nevertheless, the inherent instability of the process at high frequencies may limit its practical application.The downward continuation frequency response is defined as where h is the difference of altitude between the original survey and the target level, u and v are spatial wavenumbers in the two orthogonal horizontal directions.
The positive exponent in eq. ( 2.1) implies that an excessive amplification of high frequency noise, present in real data, cannot be avoided.Thus, in some cases the computation of the downward continued fields could not lead to meaningful results.
Several approaches have been suggested to prevent this problem, e.g.: the attenuation of the highest frequencies of the potential field by means of the smoothing process (see Dean, 1958;Grant and West, 1965); the definition of a limit value of the continuation level (Baranov, 1975); the attenuation of the amplitude spectrum only at higher wavenumbers according to a function determined by Wiener filtering theory (Pawlowski, 1995).
An alternative approach to the Fourier domain filtering is based on the approximation of the downward continued field by means of a Taylor series expansion (Evjen, 1936) (2.2) where f (x, y, z) is the potential field, z 0 is the measurement level and z the continuation level.Therefore, the computation of m-order vertical derivatives of the potential field f (x, y) is needed.To calculate any order of vertical derivatives, Fedi and Florio (2001)  Rather stable horizontal derivatives may be computed in the space domain by finite-differences or by means of splines.
Such a scheme for the computation of a downward continued field (Fedi and Florio, 2002a) involves an intrinsically stable transformation in the Fourier domain and the use of methods like finite differences for the com-putation of the second horizontal derivatives.The advantage is represented by the stability and the accuracy of the computation of the first vertical derivative, and consequently of all the higher order odd derivatives.The sum of stable vertical derivatives ensures the control of high wavenumber amplification in the downward continued field.

Case history
An example of ISVD downward continuation concerning the application of the algorithm to the real case of the gravity field of the Campanian Plain (Southern Italy) is shown here (Fedi and Florio, 2002a).The gravity data set available for the area consists of the Bouguer anomaly map (Cassano and La Torre, 1987; reduction density: 2.2 g/cm -3 ).The anomalies were digitized with a sampling step of 0.5 km (fig.2a).The Parete volcanic and geothermal area is located within the Campanian Plain (Italy).This is a Plio-Pleistocene tectonic depression filled with sediments of alluvial and volcanic origin and limited at NW, NE and SE by Mesozoic carbonate platforms and at SW by the Tyrrhenian coast (fig.1).The sinking of the carbonate platform reaches some thousands of meters (Ippolito et al., 1973).Its formation is thought to be related to the tectonic events accompanying the opening of the Tyrrhenian Sea and the anticlockwise rotation of the Italian peninsula (Scandone, 1979).These events produced a regional rise of the mantle and caused a phase of intense potassic volcanism along the Italian peri-Tyrrhenian border, specifically in some areas, as the Campanian Plain, where the carbonate basement deepens of thousands of meters.Here, volcanic complexes of Quaternary age are located (Somma-Vesuvius and Phlegrean Fields).Older evidence of volcanism was found in geothermal boreholes near Parete (calc-alkaline andesitic and basaltic lavas 1.5 km thick; Ortolani and Aprile, 1978), and northward, in the Villa Literno area.
A gravity low characterizes the Campanian area.It is justified by the sediments filling the depression.It is split into three parts by the relative highs in the Parete area and corresponding to the Somma-Vesuvius volcanic complex.The gravity high near Parete was interpreted as due to a buried volcanic structure (Baldi et al., 1976;Aprile and Ortolani, 1979) or to a buried carbonate horst with magmatic intrusions along its bordering faults, as inferred by studying the local magnetic signature (Carrara et al., 1973).The reconstruction of the buried structure generating the gravity high in the Parete area would be helped by a more detailed image of the gravity field.
The ISVD downward continuation to 0.7 km is shown in fig.2b.This continuation level was c Fig. 2a-c.a) Bouguer gravity anomaly field in the Parete area (near Naples, Italy); b) downward continuation to 0.7 km by ISVD method; c) downward continuation to 0.7 km by using the FFT linear filter.
selected after a spectral depth estimation by the method of Spector and Grant (1970) modified by Fedi et al. (1997).The maximum number of terms used in the Taylor series was five.The increase in resolution is clearly visible with respect to the original gravity field (fig.2a).For example, north of Parete two highs are now clearly singled out instead of the single positive anomaly visible in the gravity field.The northeastern one assumes an elongated shape in the NW-SE direction, while the south-western one shows now an elongation towards E and SE.The gradients between the highs and the SE low have increased considerably and are clearly W-E oriented in the southern part and N-S oriented in the northern one.Starting from these results, interpretative steps can now be greatly improved in resolution.On the contrary, the downward continuation to the same level obtained by means of the frequency response in eq.(2.1) (fig.2c) shows a large enhancement of noise in the data and, consequently, any later interpretation of this map should be necessarily preceded by a lowpass filtering step.

Multiscale analysis by enhanced horizontal
derivative method: an application to the Southern Apennines seismic active area (Italy)

Computational approaches
Potential fields are the result of the interference among anomalies of different extent presumably due to sources with different shape and extension located at different depths.Therefore, methods allowing their analysis at different scales are needed.
The Multiscale Derivative Analysis technique (MDA; Fedi, 2002) yields valuable results without resorting to any sharp separation of the effects related to frequency ranges.MDA is based on the Enhanced Horizontal Derivative (EHD) properties (Fedi and Florio, 2001).EHD may be defined as (2.3) where (2.4) and f (1) ,…, f (m) are the m-order derivatives of the field f.Finally, w 0 …w m is a set of weights.
EHD involves numerical computation of high-order vertical derivatives.To this task some stable derivation scheme is needed as the Integrated Second Vertical Derivative method (ISVD; Fedi and Florio, 2001).
The horizontal derivative of the sum of eq. ( 2.4) tends to enhance the signal over the source boundaries and to suppress the other spurious maxima due to the single derivatives of different order.Source boundaries are therefore defined by considering the location of maxima of the EHD function (eq.( 2.3)).
If a general case of sources of different depth/ extent is assumed so as to generate effects at various scales, the choice, of weight set, the starting term and the last term are decisive to obtain different images of the source boundaries.In case of correct choice, analysis of EHD for different values of the derivatives order m corresponds in practice to a Multiscale Boundary Analysis (MDA), that reinforces the potential field effects at a specific scale.MDA technique acts similarly to a filtering process, however MDA does not involve true separation of field components, but different MDA signals related to different field 'details' are determined by a selective combination of (derivative-based) components, each one having a different inherent resolution.Thus a Multiscale Boundary Analysis consists in producing a number of EHD maps having selective enhancements, and since the concept of resolution is a relative one, we will refer in the following to Large, Intermediate and Short Scale maps.

Case history
MDA should give the most interesting results where the superposition of effects due to local/ regional or shallow/deep sources, is likely to occur.Therefore, a typical test site to evaluate the efficiency of the MDA method can be re- presented by a geological domain characterized by a complex setting made-up of different structural elements, located at various depths and with local and regional extension, overlapping each other.These conditions are typically observed in systems resulting in a thick pile of tectonic units like, for example, the southern sector of the Apennine chain (Fedi et al., 2002).In fact, such a chain can be described as a complex thrust and fold belt system, built from Lower Cretaceous to Quaternary, as a consequence of the convergence between African and European plates (Finetti and Del Ben, 1986;Dewey et al., 1989;Finetti et al., 1996).
Most of the structural and geological patterns outcropping along the Southern Apennine chain and surrounding areas, which cannot be easily described by means of the simple analysis of the potential field or of other traditional methods of signal enhancement, are instead clearly re-cognized by the MDA method.However, the most significant result is represented by the existence of several structural lineaments poorly or not completely correlated to any outcropping feature, but well evidenced by MDA and with an information content often so rich as to need a very careful geological interpretation.
MDA was carried out on both gravity and aeromagnetic data sets.Gravity data windowed from the Bouguer Gravity Anomaly Map of Italy published by the CNR (Carrozzo et al., 1986; reduction density: 2.4 g/cm 3 ) were sampled with a step grid of 1 km (fig.3).Aeromagnetic field was windowed from a data set prepared by AGIP (1981), merging flight lines data to an unique altitude of around 2600 m a.s.l.(8500 feet), with a 2 km sampling step (fig.5).The example shown is relative to the area from the Tyrrhenian to the Adriatic Sea, included between 40°30′N and 41°30′N.4a), was obtained by computing EHD starting from the gravity field as first term of eq.(2.2), and considering derivatives up to m = 7 (Fedi et al., 2002).
As far as the structural setting, the tectonic style and the shallow geological features are concerned, the regional pattern of MDA reveals clearly the abrupt change in tectonic style passing from to western to the eastern side of the Central and Southern Apennines.
West of the front chain, the MDA signal shows the prevalence of short, arc-shaped and variously oriented trends, together with some linear trends with anti-Apennine direction.Most of the MDA trends within the chain seem to be related to its structural elements and generally coincide with normal faults systems and major overthrust fronts.This is consequence of the geodynamic evolution of the area, since the tectonic history of the structural units within the chain is much more complex than that at the external (eastern) side of the chain front.In fact, since Tortonian an intense crustal shortening occurred, due to a complex, poli-phasic compressive regime and to the consequent mountain building forming the current Apenninic belt.
Along the eastern side of the Apennines, several long lineaments are clearly shown with predominantly NW-SE direction.They may be related to the deformation of the external Apennines domains and should evidence a series of regional thrust fronts with linear trends.More eastward, the MDA map outlines the hidden deep normal fault systems bordering the western side of the undeformed foreland, represented by the Gargano-Apulia carbonate platform.
Gravity data: short scale map -The short scale MDA map (fig.4b) enhances the finest gravity source patterns.It was computed starting from the gravity field as first term and calculating EHD up to m = 9 (Fedi et al., 2002).The great improvement in describing structural-geologic patterns with small extension and local significance is immediately evident, although some of them do not have a clear structural meaning to be inferred on the base of the outcropping features.One of the stronger EHD maxima in this area is an uncorrelated, long regional trend extended from Abruzzo to the Gulf of Taranto, marking with high precision the position of the eastern boundary of the allochtonous chain front as inferred by well information and seismic data.
Another example is represented by a noticeable linear trend running from the Gargano to the Matera area.Its meaning could be related to the presence of an inner margin of the Apulia-Gargano Platform sunk along a regional normal fault system and buried by the Bradano Units.The southern end is partially related to structural features at the surface, where the deep carbonate block outcrops near Matera.Northward, the trend is divided in two different branches.The first one continues in the same direction and is clearly visible at the short scale.The second one turns northward and probably indicates the internal margin of a buried carbonate sector connecting the outcropping Murge and Gargano carbonate blocks.The northern segment of this branch is clearly evidenced also at the intermediate scale.
Magnetic data: short scale map -The MDA of the magnetic field of the Southern Apennines foreland area (fig.6) was obtained starting from the magnetic scalar potential as first term of the summation and calculating EHD up to m = 8 (Fedi et al., 2002).Here, the main lineaments corresponding to EHD maxima are evidenced.They are mostly NE-SW oriented (anti-Apen-ninic trend), but some of them show NW-SE trends.The lineaments that roughly follow a NE-SW trend tend to disappear just W of the front of the Appenines chain thrusts, and display an anastomosing and intersecting pattern.The absence of most of these lineaments in the gravity multiscale analysis can be explained by a small or zero density contrast along these lineaments and/or by a limited throw between the two faulted blocks.On the other hand, if faulting were ascribed to an old tectonic phase, only rocks actually lying at high depths would be affected, with a possible involvement of the crystalline basement.This may imply changes of magnetic properties, but without significant density contrasts across the fault.

Computational approaches
Euler deconvolution is one the most popular methods based on the derivative of potential fields, due to its suitability to manage several classes of sources.An 'automatic' interpretation of magnetic profiles using the Euler equation was proposed by Thompson (1982) based on early studies by Hood (1965) and Slack et al. (1967).The computation was extended to the use of map data (Reid et al., 1990) and to the gravity field and its vertical derivative (Marson and Klingélé, 1993).Some difficulties in its use in practical cases are given by the a priori selection of a constant, called structural index, that depends on the shape of the source, and by the possible existence, in a small area, of sources characterised by different structural indexes.To solve these problems, it was proposed to repeat the deconvolution several times in a moving window, varying the structural index and analysing the degree of clustering of the solutions.It was shown that in the presence of a zero background, the deconvolution problem may be solved also for the structural index (Roy et al., 2000).
The theoretical base of the Euler deconvolution method is that magnetic and gravity fields (and their horizontal or vertical derivatives) of simple sources are homogeneous (e.g., Stavrev, 1997;Blakely, 1996) (2.5) where n is the degree of homogeneity.If f is homogeneous of degree n, it may be demonstrated that .
If we consider a one-point source (that is a geometrical figure characterized by the coordinates of a unique singular point) located at (x 0 , y 0 , z 0 ) and a magnetic (or gravity) field m measured on the plane z m(x,y) = f [(xx 0 ), (yy 0 ), (zz 0 )] eq.(2.6) can be rewritten as (Thompson, 1982) f tx ty tz t f x y z In eqs.(2.7) and (2.8), N = -n represents the structural index.Thompson (1982) pointed out that the presence of regional fields or dc offsets masking the absolute level of an anomaly, makes unstable the application to real data.Thus Thompson (1982) rearranged eq. ( 2.8) considering that the measured potential field is the sum of the anomaly field and a background: (2.9) The knowledge of the anomalous field and of its derivatives in eq.(2.9) leads to a linear system of equations.In the most common approach (Thompson, 1982;Reid et al., 1990), this system may be solved in a moving overlapping window for the unknown coordinates of the anomaly source (x 0 , y 0 , z 0 ) and B, assuming an a priori knowledge of the structural index N.The inversion is repeated using different values for N until the tightest clustering of solutions for a particular feature is finally selected.Serious drawbacks must however be mentioned: -Existence of sources with shapes intermediate between two simple models, thus having fractional structural indexes with a value variable with the distance between source and observations (e.g., Slack et al., 1967;Steenland, 1968).
-Implicit assumption that in each analysed window only the effect of a single source is present.
-Variation of N depending on the analysed single portion of the considered anomaly and, consequently, influence of the direction of the total magnetization vector (in the magnetic case) on the parameters estimation (Ravat, 1996).
-Impossibility to simultaneously estimate the depth source and N, because they are linearly dependent.
A new approach to Euler deconvolution was recently suggested by Fedi and Florio (2002b).The technique consists in the application of eq. ( 2.8) to a n-order vertical derivative of the field, instead of the field itself (Derivative Euler Deconvolution, DED).The application of eq.(2.8) to a n-order vertical derivative suppresses the constant background level, reduces long period trends and, increasing field spatial resolution, reduces the problems of interference between nearby anomalies.As a consequence, this allows eq.(2.8) to be solved simultaneously both for the position of the sources and the structural index, giving additional information about geologic structures.In this way, use of Euler equation is generalized to any source model and there is no need to assume any a priori (and then subjective) value for the structural index.Such an approach can be very useful in complex structural regions.
Obviously, the application of this method also implies the computation of k + 1 order horizontal and vertical stable derivatives.The problem can be solved by means of stable algorithms such as the Integrated Second Vertical Derivative method (ISVD, see Section 2.1; Fedi and Florio, 2001).

Case history
The Derivative Euler Deconvolution provides a classification of the sources basing on their shape, so giving valuable information for the interpretation of complex geologic structures like a volcanic district.A volcanic environment is in fact characterized by very different structure typologies causing gravity effects (faults zones, lithology contacts, igneous bodies, dikes, etc.).Phlegraean Fields (fig. 1) is an area of active volcanism west of the town of Naples (Italy).In this area, the volcanic activity of the last 37 000 years was mainly explosive.Most of the products are pyroclastic rocks, with composition ranging from trachybasalts to phonolitic alkali-trachytes (e.g., Rosi and Sbrana, 1987).Great thickness of light pyroclastites was confirmed (San Vito wells) in the central part of the district by geothermal exploratory well logs (Cassano and La Torre, 1987), together with the lack of thick lava complexes that, instead, are present along the peripheral sectors (AGIP, 1987).This supports the hypothesis of a caldera formed as a consequence of major eruptions which occurred in this area (37000 years b.p.Campanian Ignimbrite eruption; 12 500 years b.p., Neapolitan Yellow Tuff eruption).Since about 5 ka, besides a recent minor eruption (Monte Nuovo eruption, 1538), the major tectonic event is represented by a series of uplifts inside the caldera.The last ones occurred in 1970-1972 and 1982-1983, centered in the area of the town of Pozzuoli (i.e.Orsi et al., 1996) with a net rise of about 3.5 m.These uplifts were interpreted as the consequence of the arrival of new magma into a shallow reservoir (Bonafede, 1990;De Natale et al., 1991) or as the combined effect of the latter with a pore pressure increase (Casertano et al., 1976).
Gravity data were digitised from the AGIP (1987) Bouguer anomaly map data (reduction density: 1.4 g/cm 3 ; Barberi et al., 1991).The most significant feature is represented by a circular gravity low centred on the Gulf of Pozzuoli and surrounded by a belt of positive anomalies.An intense positive anomaly is present in the SW area between Procida and Monte di Procida, while the lower values are located in the NE area.The existence of a main structure reflecting the geometry of a collapsed caldera structure, was already inferred by Florio et al. (1999) by analysing the horizontal derivative of the gravity and magnetic maps of Phlegraean Fields.Other structures are present inside the caldera, and a good correspondence with epicenters of recent seismicity disclosed young ages for these features.Other lineaments probably having regional importance were also found, with SW-NE direction.
DED was applied (Fedi and Florio, 2002b) to the vertical derivative of gravity field (fig.7a) to recover 3D coordinates of sources (fig.7b) and structural indexes (fig.7c).DED clearly delineates the SW sector of the inner boundary of the caldera.In the western sector it shows structural indexes typical of limited-throw faults (0.4 < N < 0.7), the top of which is computed to a depth of about 0.15 km.In the eastern sector, the structural index has higher values and the top of this structure seems to be at about -0.2 km of depth.The same model (limited-throw faults) is confirmed for the most part of the outer boundary of the pericalderic magmatic intrusions, where computed depths range between -0.1 and -0.35 km.In correspondence with the belt of gravity highs, a number of structures having a higher N value (sheet-like sources, 0.75 N < 1.75) show up, consistently with the fact that in this area sources are expected to be sills.Solutions associated with the highest structural index values (linear or point like sources, 1.75 N < 3.0) concentrate in a few points of the map that correspond to known volcanic vents.This is consistent with the probable presence of dense igneous bodies along the conduit.

Boundary analysis of the 3D gravity
anomaly field of Southern Apennines active seismic area

Description of the method
As broadly treated in the previous sections, the study of derivatives of the potential fields represents a signal enhancement technique experimented since a long time.It allows the information content inside the signal to be amplified without implying any arbitrary assumption.
The simplest procedure is represented by the analysis of the horizontal derivative of the potential field, defined as follows: (2.10) where M is the potential field and h is its horizontal derivative.
The horizontal derivative assumes the highest values just in correspondence with the lateral boundaries of the source anomaly.This additional information, not immediately evident from the natural field, helps to provide a reconstruction of the geological setting of the structures investigated.However, this technique is fully reliable only in the case of vertical or sub-vertical discontinuities.In fact, in the presence of oblique surfaces, maxima of the horizontal derivative are located on the vertical correspondence of a point of the discontinuity located at an intermediate depth (Grauch and Cordell, 1987), without any information on the dipping direction.
Therefore, the efficacy of this simple method is good in the case of structural features typical of an extensional regime (normal faults, transform faults) or in the presence of dykes, volcanic conduits, plutonites, etc.On the contrary, it is scarcely exploitable when the structural setting is characterised by a compressive tectonic style (presence of inverse faults, transpressive faults, multiple thrust systems, etc.).
A strategy to improve this method, consisting in several steps, was proposed (Cella et al., 2000): a) application of a Hanning tapering window to the potential field to prevent edge effects in the following step; b) upward continuation to different increasing levels of the potential field; c) computation of the horizontal derivative; d) automatic location and thresholding of maxima values of the gradient of the field continued at different levels, following the method suggested by Blakely and Simpson (1986).
If the anomaly field, due to a density contrast along an oblique discontinuity, is upward continued at increasing levels, the resulting signal will be related to progressively deeper sectors of the discontinuity.Therefore, the position of the maxima of the horizontal derivative of the continued field will be laterally shifted toward the dipping direction of the discontinuity, proportionally to the continuation level.
Synthetic tests, made on several kinds of simple sources, confirmed this prediction and showed that the shifting rate of maxima of the gradient is subjected to an asymptotic decrease with increasing continuation levels.Thus, major shift is observed at lower continuation levels.Only in the presence of a regional trend or in complex cases (when several significant sources lay very close), the interference can be so intense as to hide the described effect.Removing the long period component of the signal or, in the second case, decreasing the step grid of the data set can help to obtain good results.This allows the range of upward continuation levels to be concentrated to the lowest values, where the effect due to the dipping sense prevails.
Another strategy is the use of the vertical derivative instead of the horizontal one.The vertical gradient in fact allows a further increase in resolution of the local anomalies due to a higher signal decay with distance.
Another positive feature of 3D boundary analysis is that it allows a fast qualitative evaluation of the depth to sources basing on the progressive disappearance of the trends of maxima of the horizontal gradient increasing the continuation level.The higher the level at which the trend disappears, the deeper the source.

Case history
A preliminary application of the method on the Bouguer gravity field of the Central-Southern Apennine is presented (Cella et al., 2000).To prevent the effect due to the presence of a regional field, data were filtered by using a high-pass numerical filter.The residual field (fig.8) at several continuation levels was horizontally derived and maxima were plotted as points with different colours depending on the continuation level (fig.9).
As predicted by several tests, the results consist in many trends of maxima variously shaped and showing a different shifting with increasing continuation levels.This allows vertical discontinuities (no shifting) to be distinguished from oblique surfaces and, in this case, to recognise their dipping direction.Consequently, many complex structural settings can be identified along the Apennine belt (Cella et al., 2000).

Conclusions
The spatial derivative-based techniques discussed above are included in a large number of methods that do not imply a data interpretation by means of classical approaches (forward or inverse method), but enhance the gravity and magnetic signals to extract a large amount of information inherent in such data sets, but not directly evidenced by the measured fields.This information generally concerns the horizontal position of the boundaries separating the anomaly sources.As the process of derivation is inherently unstable, to compute high order derivatives, it is necessary to use specially designed algorithms, such as the Integrated Second Vertical Derivative (ISVD; Fedi and Florio, 2001).
In this paper we illustrated several improvements of classical enhancement techniques.
These improvements regarded: 1) A stable downward continuation scheme, based on a Taylor series approach and on the computation of ISVD vertical derivatives.
2) A multiscale boundary analysis (MDA), based on the Enhanced Horizontal Derivative technique.By MDA it is possible to produce a number of boundary maps related to different scales without any true separation of field components.
3) A version of the Euler Deconvolution (Derivative Euler Deconvolution) method in which there is no need to assume a model of source, but instead the type of source is one of the searched parameters, together with the 3D position of sources.
4) A 3D version of the horizontal derivative method, where the boundary analysis is performed on a set of upward continued fields.This can give additional information, specifically about the inclination of density/magnetization contrasts.
(ISVD) method, by which the first vertical derivative of the potential field is computed in two different steps: a) vertically integrating the field by means of a frequency domain operator; b) computing the second vertical derivative of the integrated field by means of the sum of its second horizontal derivatives, according to the Laplace equation.Other vertical derivatives are computed by the Laplace equation starting from the field (any order even derivatives) and its first vertical derivative (any order odd derivatives).

Fig
Fig. 4a,b.a) Intermediate scale MDA map of Bouguer anomaly gravity field in the area windowed as shown in fig. 3. b) Short scale MDA map of Bouguer anomaly gravity field in the area windowed as shown in fig. 3. Lineaments are represented by trends of maxima, which are graphically recognised by the colorbar in normalized units.

Fig. 6 .
Fig. 6.Short scale MDA map of magnetic data field in the foreland area windowed as shown in fig. 5. Lineaments are represented by trends of maxima, which are graphically recognized by the shown colorbar in normalized units.

FigFig. 8 .
Fig. 7b,c.b) DED solutions shown as function of their depth; c) DED solutions shown as function of the related structural index.The found depths, associated with N estimates between 0 < N < 3, are in the range + 0.12 > z 0 > -1.28 km.Gray shading represents topography (darker colours = higher values).Coordinates UTM, zone 33N.

Fig. 9 .
Fig. 9. Trends of local maxima of the horizontal gradient of the field shown in fig.8, upward continued to several increasing levels (from 1 to 30 km).