GEOPHYSICS, VOL. 71, NO. 1 (JANUARY-FEBRUARY 2006); P. L5–L11, 14 FIGS.

10.1190/1.2163912

Leveling aerogeophysical data using a moving differential median filter

Eirik Mauring1 and Ola Kihle1

Their crossline wavelength is often twice the flight-line spacing

ABSTRACT (Minty, 1991). Most microleveling or decorrugation routines
involve filtering a grid of data that is orthogonal to the line
We describe a new technique that can be used to level direction to detect remaining errors along lines (Luyendyk,
data collected along regular and irregular line patterns 1997). A common approach used by major software vendors
with or without tie-line control. The technique incor- combines a 2D high-pass filter and a directional filter to mi-
porates a moving differential median filter to minimize crolevel data (Ferraccioli et al., 1998). A similar microleveling
line-level errors, to level survey-line data, and to micro- technique for gridded data is described by Urquhart (1989).
level data with no tie-line control. This overcomes the The disadvantages of this approach are that resolution is lost
problem of standard leveling methods that lose their ef- during the preliminary gridding step, and the approach is not
fectiveness with irregular flight patterns. designed for survey data with irregular lines and line spacings.
To validate the method, we use it to level very-low- This paper introduces a microleveling/leveling technique
frequency (VLF) electromagnetic (EM) data from a he- that can be applied to data collected along both irregular and
licopter survey where flight lines are parallel. Leveling regular line patterns with or without tie-line control. The tech-
is also performed on a set of vintage aeromagnetic data nique is a modification of a technique suggested by Liukkonen
from the North Sea, gathered from nonparallel flight (1996), who uses it to process radiometric data. Both direc-
lines. Results show that the differential median filter tional and differential median filtering can be used for leveling
leveling technique is superior to the standard leveling and microleveling aerogeophysical data. Microleveling is nor-
method because it results in fewer line errors and less mally the term used for small line-error corrections applied
distortion of high-wavenumber anomalies when proce- after tie-line leveling (Luyendyk, 1997). None of the examples
ssing irregular survey lines, making the method suitable in this paper have tie-line control, so leveling is the most ap-
for a wide variety of data sets. propriate term to use for the differential median filter leveling
technique.

INTRODUCTION LEVELING TECHNIQUES

Aerogeophysical data often suffer from line-level errors This section compares a standard leveling technique with
(corrugations), which are data inconsistencies between ad- the differential median filter leveling technique. Our proposed
jacent flight lines arising from a variety of unavoidable technique is shown to be superior to standard leveling.
acquisition-related circumstances. This problem is partly
solved by flying in a grid pattern with lines and tie lines and Standard leveling technique (directional filtering)
then resolving the inconsistencies during data processing.
However, data collected during difficult conditions [e.g., time- The directional filtering technique uses sampled values
varying fluctuations in the earth’s magnetic field, rough ter- from an error grid that (ideally) shows line errors attributable
rain, unstable very-low-frequency (VLF) transmitters] may to insufficient leveling. A grid data set is generated that best
still show corrugations after various corrections and tie-line describes line-level errors in survey data. Traditionally, this
leveling. is done by applying a high-pass filter and a directional fil-
The last step in the leveling procedure may be to apply ter to gridded data that have first been subjected to various
a microleveling routine to remove the remaining line-level corrections and tie-line leveling (Luyendyk, 1997; Ferraccioli
errors. Such routines are based on the assumption that line er- et al., 1998). It is customary to use a high-pass filter length
rors appear as elongated artifacts along flight lines (Urquhart, of four times the line spacing on the preprocessed grid and
1989; Minty, 1991; Luyendyk, 1997; Ferraccioli et al., 1998). subsequently to apply a carefully designed directional filter to

Manuscript received by the Editor January 31, 2003; revised manuscript received March 22, 2005; published online December 19, 2005.
1
Geological Survey of Norway, 7491 Trondheim, Norway. E-mail: [Link]@[Link].


c 2006 Society of Exploration Geophysicists. All rights reserved.
L5
L6 Mauring and Kihle

ensure that only energy nearly orthogonal to the lines is The technique is based on filtering line data, as opposed to fil-
passed. Using this combination of high-pass and directional tering grid data when using directional filtering. For a given
filters is known as directional filtering. line, the 1D median is determined at each station based on
The combination of the high-pass and directional filters en- data values along the line within a given distance d/2 from the
hances short-wavelength features across the line direction, station (see Figure 1). In a similar way, the corresponding 2D
which are considered to be line errors. The resulting grid is median value is determined from nearby data values in the
sampled along the survey lines. The sampled values are then current line and inside the circle that intersects neighboring
smoothed using a nonlinear filter followed by low-pass filter lines. The difference between the 2D and 1D median values is
convolution (Fraser et al., 1966) to avoid removing any re- taken to be the leveling error, and the level correction is made
maining geologic signal. These correction values are then sub- by adding this value to the data value at the current station.
tracted from the preprocessed values to produce the final lev- With the directional filtering technique, it is important that
eled data. the corrections are smoothed to avoid removing any remain-
ing geological signal (see Figure 3).
For stable results with median filtering, a regional trend
Differential median filter
must first be removed from some types of data (Smith and
Figure 1 shows 1D and 2D windows for different line pat- Wessel, 1990). Two standard methods are 2D high-pass filter-
terns used in median filtering. A circular window with radius ing and trend removal based on a polynomial function. After
r is used to determine the 2D median. A corresponding 1D
filter with window length d is placed over the central line seg-
ment in the middle of the 2D median window. The point to be
processed is centered within the line segment. The windows
pass over the data set in unison, station by station, generating
level-correction values, one at a time.
For the 1D filter, the median is defined as the midpoint in
an array of numbers sorted in ascending order. The median is
denoted

median (X1 , X2 , . . . , Xn ).
The 1D median filter of size n (odd) in the sequence {Xi ,
i ∈ Z} is defined as (Justusson, 1981)


yi = median xi = median (xi−v , . . . , xi , . . . , xi+v ), i ∈ Z,
n

where v = (n − 1)/2 and Z denotes all natural numbers. The

nonlinear median filter has two very important properties: It
acts as a noise filter, and it is particularly efficient in remov-
ing spikes. Although it is a smoothing filter, it preserves sharp
edges in a data set (Gallagher and Wise, 1981; Justusson, 1981;
Stewart, 1985). Both properties are illustrated in Figure 2. The
figure also outlines the 1D median filtering procedure.
A 2D median filter is a simple extension of the 1D filter to
Figure 2. Procedure for the median filtering of lines using
an array of numbers contained within an area. A 2D median a five-point filter (after Stewart, 1985). Note the removal of
filter with window A in a data set {Xij (i, j ) ∈ Z 2 } is defined by spikes and the preservation of edges in the data.


yi = median xij = Median [xi+r,j +s ; (r, s) ∈ A],
A
i, j ∈ Z 2 . (1)

Figure 1. The 1D and 2D windows used in differential median

filtering for different line patterns. The length of the 1D filter Figure 3. Smoothing raw correction values using a Naudy filter
is d, and the radius of the circular 2D window is r. preceding the application of a low-pass filter.
 Moving Differential Median Filter L7

determining correction values, we add the smoothed correc-

tion values to the preprocessed anomaly field line data. These
data are then regridded.

TESTING LEVELING TECHNIQUES ON

AEROGEOPHYSICAL DATA
In this section, two survey examples are presented to eval-
uate the quality of the differential median filter in compari-
son to directional filtering. Neither survey includes tie lines.
The two examples demonstrate cases in which (1) lines are
parallel but the area is irregular and (2) lines are irregular in
spacing and direction. The parallel-lines example uses VLF
EM data from a helicopter survey. The example of irregularly
flown lines uses vintage offshore aeromagnetic data acquired
by the Geological Survey of Norway over the North Sea in the
1970s (Figure 4).

Helicopter survey
A map of the flight lines is shown in Figure 5. The purpose
of the survey was to provide information for mineral explo-
ration. The methods used were electromagnetic, VLF EM, ra-
diometric, and magnetic. The line spacing was 200 m. The fly-
Figure 4. Overview map with flight lines for the vintage aero-
ing height was 60 m, and lines were flown in south and north magnetic survey.
directions (Figure 5). A minimum curvature routine (Briggs,
1974) using a grid-cell size of 50 m was applied when gridding
the data.
We have chosen the results from leveling the VLF EM data.
The VLF instrument receives the magnetic component of
fields radiated from one or two VLF transmitters (used
around the world for navigation and to communicate with sub-
marines). The change in total field (real component) is mea-
sured. The instrument has a sensor comprised of three mutu-
ally orthogonal antennas that record signals from transmitters
radiating from different directions. This provides good cou-
pling with ground conductors of any orientation (RMS Instru-
ments, 2002). In this paper we show the real-component re-
sults from the inline antenna that couples maximally with a
field in the nominal direction of flight (inline with flight). This Figure 5. Map of the survey area with flight lines for the he-
field was transmitted from a 16-kHz VLF station transmitter licopter survey with the direction of the traditional leveling
in Rugby, England. filter.
The power output could not be predicted or
controlled during the survey. Therefore, a con-
stant base level was estimated and removed,
leaving lines containing modified anomaly
magnitudes. The grid-cell size was 50 m, and
the sampling interval was approximately 5 m.
Gridded, base-removed data are shown in
Figure 6a. Line-level errors are evident. These
are mainly caused by fluctuations in the VLF
EM data from unaccounted elevation changes
in the flight patterns. The percentage values
on the legend bar are deviations from a nomi-
nal 100% primary magnitude determined in an
anomaly-free zone (RMS Instruments, 2002).
The data were leveled using both directional
microleveling and differential median level- Figure 6. (a) Grid of helicopter survey data after preprocessing but before lev-
ing. Directional leveling was performed using eling. (b) Line-error grid after 800-m Butterworth high-pass and 0◦ directional
a pass angle of 0◦ . The line-error grid is shown filtering.
L8 Mauring and Kihle

in Figure 6b. An 800-m Butterworth high-pass filter was ap-

plied to the grid to remove trends in the differential median
leveling. Here, we used a 1D filter length of 200 m and a 2D
filter radius of 450 m. We tried 2D filter radii of 250 m (1.5
times the maximum line spacing), 450 m (2.5 times the maxi-
mum line spacing), 650 m (3.5 times the maximum line spac-
ing), and 850 m (4.5 times the maximum line spacing). The
absolute difference between the correction values for the dif-
ferent radii were 0.51, 0.21, and 0.13 nT, respectively, while
the 1D filter length was kept constant (200 m). Figure 7 shows
these differences as curves for one of the survey lines. The dif-
ferences between the 450 to 650-m and the 650 to 850-m radii
are negligible. The absolute differences indicate that using a
2D window larger than 650 m (3.5 times the line spacing) has
little effect on the quality of the final leveled result.
Comparing the two grid images in Figure 8 reveals practi-
cally no line errors, neither for the traditional leveled data nor
Figure 7. Curves showing differences in correction values for for the differential median leveled data. To highlight possi-
varying 2D filter radii. The 1D filter length is kept constant ble line errors, the first vertical derivative of the anomalies is
(200 m). shown in Figure 9. Presenting the first vertical derivative of the
VLF magnetic component can be justi-
fied because the magnetic field is a poten-
tial in the air in the quasistatic approxi-
mation (Pedersen and Oskooi, 2004).
The grids of the first vertical derivative
show no remaining line errors, only high-
wavenumber anomalies from electrical
conductors (geologic structures or man-
made installations such as power lines and
fences). Comparing anomaly profile curves
did not reveal any distortions of anoma-
lies either. Thus, there appear to be no no-
ticeable differences in the performance of
the two leveling techniques for this survey
area.

Vintage offshore magnetic survey

Figure 8. (a) Grid of directionally leveled data. (b) Grid of differential median lev-
eled data with a 1D filter length of 200 m and a 2D filter radius of 650 m. The correc- A detailed map of the survey area with
tion values for both techniques were smoothed using a 100-point Naudy filter and a flight lines is shown in Figure 10. No tie
1-km low-pass convolution filter applied to the preprocessed data. lines were flown, but adjustment was made
for base magnetometer data, despite the
fact that the area is situated too remotely
from the base magnetometer station for
the station readings to be used to reduce
of line errors. According to Reeves (1993),
the most remote part of the survey area
should be within 50 km distance from the
base magnetometer station to be used in
this way. The fan-shaped line pattern re-
sults from the use of Decca lanes for navi-
gation. The survey was carried out using a
maximum line spacing of 8 km. The sam-
pling interval along the lines is 150 m and
the flight altitude is approximately 500 m.
In addition to base magnetometer cor-
rections, preprocessing included spike re-
moval and IGRF (1975) corrections. Grid-
Figure 9. (a) First vertical derivative of the VLF anomalies of (a) directionally lev- ded preprocessed data are shown in Fig-
eled data and (b) differential median leveled data. On both grids there appear to be ure 11a; line-level errors (corrugations) are
no remaining level errors. evident.
 Moving Differential Median Filter L9

The data were leveled using both directional leveling and error adjustments to work properly. Some guidelines for set-
differential median leveling. First, a 35-km Butterworth high- ting these parameters are outlined below.
pass filter was applied to the grid. This is equivalent to about The length of the 1D filter along the lines (distance d in Fig-
four times the maximum line spacing. This high-pass filtered ure 1) should be less than twice the length of any level er-
grid is the basis for both leveling techniques. For the median rors along the line, according to Liukkonen (1996). However,
filtering this was done for trend removal. Directional leveling removing all level errors can distort anomalies (shape and
was performed using a pass direction at an angle of 30◦ clock-
wise from the south. This is the average direction of the survey
lines. The resulting line-error grid for directional leveling is
shown in Figure 11b. For differential median leveling, the 1D
filter length is 20 km and the 2D filter radius is 28 km (3.5 times
the maximum line spacing). We tried radii of 12, 20, 28, and 34
km while keeping the 1D filter length constant. The absolute
average differences between the correction values were 4.4,
1.5, and 0.9 nT, respectively, for all of the flight lines. As for
the helicopter survey, using a filter radius larger than 3.5 times
the maximum line spacing seemed to have negligible effect on
the quality of the final, leveled result. Further processing steps
were carried out as depicted in the earlier sections discussing
leveling techniques.
Comparing the two grid images in Figure 12 reveals obvi-
ous line-level errors, especially in the northwestern part of
the directionally filtered grid (Figure 12a). This is best seen
in the corresponding high-pass filtered grid image of Figure 13
(35-km Butterworth high-pass filter). This is from the large an-
gle deviation from the 30◦ directional filter. Also, some high-
wavenumber anomalies are attenuated on the directionally
leveled grid (see circled areas in Figure 12a). However, small
corrugations can also be seen in the southeastern part of the
differential median leveled grid (Figure 13b), although less
prominent than for the directionally leveled grid. Figure 10. Map of the survey area with flight lines for the vin-
Figure 14 shows magnetic anomaly values for lines P4 and tage aeromagnetic survey with the direction of the traditional
P58 (see Figure 10 for location) before and after leveling with leveling filter and the positions of lines in Figure 14.
the two techniques. The figure reveals se-
vere distortion of high-wavenumber anoma-
lies for the directional leveling technique in
areas where the line direction deviates from
the optimal direction for the filter. The dis-
tortions are most prominent between 55 and
70 km on line P4. On line P58, the distor-
tions are most severe between 40 and 60 km
and between 120 and 140 km. It is evident
that the differential median filter is superior
to the directional filter in leveling the vin-
tage data with less distortion of anomalies
and independence of line direction.

DISCUSSION
This section discusses how to set 1D and
2D parameters. It also gives examples of
computing time for the case histories pre-
sented.

Parameter settings
A disadvantage of the differential median Figure 11. (a) Grid of vintage data after preprocessing (spike removal, base mag-
filter is that the parameter settings (espe- netometer, and IGRF corrections) before leveling is applied. (b) Line-error grid of
cially the 1D filter) need some trial-and- vintage data after 35-km Butterworth high-pass and 30◦ pass angle filtering.
L10 Mauring and Kihle

amplitude). Subsequent modeling of these

anomalies (e.g., calculations of depth to
magnetic basement) may then become er-
roneous. The gridded data may, however,
look smooth when displayed on a map. It is
crucial that we find the optimal parameters
for the 1D and 2D filters. In addition, the
optimal parameters may vary within the
survey area. The challenge is to find the pa-
rameters that preserve the shape and am-
plitudes of the anomalies while giving a
smooth appearance to the gridded data. A
trial-and-error approach can be done for
the 1D filter length.
For the 2D filter, we suggest that its ra-
dius (r in Figure 1) should be larger than
the maximum distance between at least
three neighboring lines to each side of the
line being processed. For the two survey
examples presented, the corrections con-
verge toward similar values as the 2D ra-
dius increases beyond three times the max-
imum line spacing. This is a rule of thumb
derived by comparing the leveling correc-
Figure 12. (a) Grid of directionally leveled data. (b) Grid of differential median tion values for different sizes of the 2D
leveled data with a 1D filter length of 20 km and a 2D filter radius of 28 km. For window while the length of the 1D window
both techniques the level correction values are smoothed using a 20-point Naudy
filter and a 10-km low-pass convolution filter applied to the preprocessed data. Note is kept constant. It is imperative that the
attenuation of high-wavenumber anomalies on (a) (circled areas). corrections be smoothed to avoid remov-
ing any remaining geological signal (see
Figure 3).
The variety of parameter settings gives
more flexibility for skilled processing per-
sonnel than the standard technique, which
is quite straightforward and does not
allow for much variation of processing
parameters.

Computational time
The new technique requires longer pro-
cessing time on a computer than the tra-
ditional microleveling method, and it in-
volves sorting a 1D and 2D array for each
data point to be processed. Computing
time is dependent on the number of data
points, data-sampling intervals, and size of
the 1D and 2D filters. However, with ever-
faster computers developing, this should
not be an obstacle for using the technique.
To filter the data for the helicopter survey,
it took 7 minutes and 45 seconds to process
315 914 points with a sampling interval of
5 m using a 1D filter length of 200 m and
a 2D filter radius of 650 m. For the vintage
offshore survey, it took 1 minute to process
Figure 13. (a) Residual grid of directionally leveled data. (b) Residual grid of differen- 120 740 points with a sampling interval of
tial median leveled data. There are fewer remaining corrugations on the differential 150 m using a 1D filter length of 20 km and
median leveled data. Both grids are produced after applying a 35-km Butterworth a 2D filter radius of 28 km.
high-pass filter.
 Moving Differential Median Filter L11

helicopter survey with a regular line pattern

and vintage aeromagnetic data collected in
a fan-shaped irregular line pattern. For the
vintage survey, the new technique effec-
tively removes line-level errors with less
anomaly distortion than the directional fil-
tering technique. We suggest that the new
technique can be used for many types of
data.

ACKNOWLEDGMENTS
The authors thank Mark A. Smethurst
for valuable comments and corrections to
this paper. They are grateful to Crew De-
velopment Corporation for data from the
helicopter survey. Thanks also go to V. J. S
Grauch and R. O. Hansen for reviewing
this paper and providing suggestions for im-
proving it.
Figure 14. Magnetic anomaly values before and after leveling lines P4 and P58 of the
vintage data. Distortion of high-wavenumber anomalies on the directional filtered
data are evident.

REFERENCES
Briggs, I. C., 1974, Machine contouring using minimum curvature:
General remarks Geophysics, 39, 39–48.
Ferraccioli, F., M. Gambetta, and E. Bozzo, 1998, Microlevelling pro-
We have applied the technique with success on a wide vari- cedures applied to regional aeromagnetic data: An example from
ety of data, from small helicopter surveys to large gravity and the Transantarctic Mountains: Geophysical Prospecting, 46, 177–
196.
magnetic compilations. While the directional filter works on Fraser, D. C., B. D. Fuller, and S. H. Ward, 1966, Some numerical
regular, parallel lines, the differential median filter seems to techniques for application in mining exploration: Geophysics, 31,
1066–1077.
work also on irregular and nonparallel lines with less distor- Gallagher, N. C., and G. L. Wise, 1981, A theoretical analysis of
tion of anomalies. the properties of median filters: IEEE Transactions in Acoustics,
Speech and Signal Processing, 29, 1136–1141.
Justusson, B. I., 1981, Median filtering: Statistical properties, in T. S.
Huang, ed., Two-dimensional digital signal processing I. Trans-
forms and median filters: Springer-Verlag.
Liukkonen, J., 1996, Levelling methods for aerogeophysical data:
CONCLUSIONS 10 August 2005. [Link]
html.
Two leveling techniques applied to remove line errors (cor- Luyendyk, A. P. J., 1997, Processing of airborne magnetic data: Jour-
rugations) from a set of aerogeophysical data are compared. nal Australian Geology and Geophysics, 17, no. 2, 31–38.
These are based on directional filtering and on a new tech- Minty, B. R. S., 1991, Simple micro-levelling for aeromagnetic data:
Exploration Geophysics, 22, 591–592.
nique built on differential median filtering. The first technique Pedersen, L. B., and B. Oskooi, 2004, Airborne VLF measurements
is widely used in the geophysical data-processing industry. It and variations of ground conductivity: A tutorial: Surveys in Geo-
physics, 25, 151–181.
involves the combined use of 2D high-pass and directional fil- Reeves, C. V., 1993, Limitations imposed by geomagnetic variations
ters, and it requires that lines be parallel. An alternative lev- on high quality aeromagnetic surveys: Exploration Geophysics 24,
eling technique is suggested that utilizes a moving differential 115–116.
RMS Instruments, 2002, HERZ TOTEM-2A VLF electromagnetic
median filter. As opposed to traditional leveling techniques, system user’s guide: GRMS Instruments.
the alternative technique can be used for nonparallel flight- Smith, W. H. F., and P. Wessel, 1990, Gridding with continuous cur-
lines and irregular flight-line patterns. For some data types, vature splines in tension: Geophysics, 55, 293–305.
Stewart, R. R., 1985, Median filtering: Review and a new F/K analogue
the technique requires removing a regional field before it can design: Journal of Canadian Society Exploration Geophysicists, 21,
be applied. 54–63.
Urquhart, T., 1989, Decorrugation of enhanced magnetic field maps:
We have tested the new technique on a wide variety of data 59th Annual International Meeting, SEG, Expanded Abstracts,
sets with success, two of which are presented as examples: a 371–372.