Let us consider a simple 2D geological model (x–z plane) composed of a diffraction point (diffractor) placed on a perfectly resistive medium with a constant EM velocity. The coordinates of this diffractor are x_d and z_d, respectively (Fig. 1a). We assume a zero-offset survey with transmitting and receiving antennas that move on a flat horizontal surface at z = 0 (dashed line in Fig. 1a). In this case, the result of the zero-offset GPR profile in time (x–t plane) will be a diffraction hyperbola (shown by the dashed line in Fig. 1b) and the electric field variation can be described by a scalar wave propagation equation, which is similar to the acoustic wave equation (Leparoux et al., 2001).

The goal of the migration is to find the geological model (in the x–z plane) from the zero-offset GPR profile (in the x–t plane). For a resistive medium (high-frequency approximation), we can use Kirchhoff's method that gives the wave field at the location of the diffractor (x_d, z_d) from the zero-offset wave field measured at the surface z = 0 (Feng et al., 2009; Schneider, 1978). Practically, the Kirchhoff migration will calculate the diffraction hyperbola (migration template) for each point of the GPR profile and, by adding the amplitudes along the template, will place it at the top of the template in the migrated profile (Claerbout, 1985; Yilmaz, 2001). Migrating each of these points for a given velocity will focus the amplitudes at their correct positions and the reflector is imaged with its true position and dip angle.

When GPR measurements are performed over a surface with topography, the migration template is no longer a diffraction hyperbola; instead, it will be a distorted diffraction curve. This is shown in Fig. 1 where the topography is chosen to be a circle whose centre is on the diffraction point (x_d, z_d), and on both sides of which the topography is flat (see the thick line in Fig. 1a). The distance between the diffraction point and the antennas (in zero offset), moving on the surface along the circle, is constant. Therefore, the migration template, shown by the thick line in Fig. 1b, will be flat on the top, and on both sides it will be represented by two flanks of a diffraction hyperbola. In this case, the imaging result of the classical Kirchhoff migration with a flat datum plane will be spurious (Fig. 4). For this reason, we absolutely need to take into account the topography of the GPR acquisition surface.

For the standard Kirchhoff migration, at a location x on the surface z = 0 (i.e. the antennas move on the flat datum plane), the two-way travel time t(x) along the grey line path in Fig. 2 is given by:

Correcting for the topography means to choose for the migration template the thick line of Fig. 1b, instead of using the dashed one, which is exactly a diffraction hyperbola. This will allow the template to follow exactly the real travel path of the GPR data. Indeed, for the same x location (Fig. 2), the z position of the antennas (moving on the rugged surface) has been changed and the two-way travel time t(x) is now calculated along the thick line path in Fig. 2 to obtain:

with:

This equation is the same as the one given without any demonstration by Lehmann et al. (1998) in the case of a 2D migration. Fig. 3 presents an overview of the different steps to compute the topographic migration. We have used the same notations as in equations (1) to (4).

In Fig. 4 we compare the results of the classical migration with flat datum plane, the classical migration after static shift and the topographic Kirchhoff migration. Fig. 4b displays the synthetic radargram computed with the model of Fig. 4a. This radargram is obtained by using a second order Ricker source having a dominant frequency of 500 MHz, located over a homogeneous medium with a velocity of 0.1 m/ns. The distance between traces is 2 cm.

Fig. 4c shows the classical migration of the zero-offset GPR synthetic data of Fig. 4b. One can see a flat horizontal 2-m-wide layer located at a depth of 1.5 m, as well as a bright spot in the middle of the section at a depth of around 1 m (Fig. 4c). The imaging result is very poor and might lead to a misinterpretation of the data. The actual classical procedure is a static shift followed by a classical migration. Fig. 4d shows the synthetic data after the static shift, and Fig. 4e displays the migration after the static shift. The result seems to be better than the one of Fig. 4c. In Fig. 4e we observe not only a bright spot at the correct depth of 1.5 m, but also two strong spots located on both sides of the diffraction point (around a depth of 1.2 m). In Fig. 4f, we present the result of the topographic migration appropriately weighted by an amplitude factor proportional to cos(α) = t_top/t(x) which also depends on the topography (Fig. 2). The amplitude factor is added to take into account the directivity factor that describes the angle dependence of amplitudes and is given by the cosine of the angle between the propagation direction and the vertical axis (Claerbout, 1985; Yilmaz, 2001). The data are well imaged and, as expected, are focused on a single bright spot located at its real depth of 1.5 m (Fig. 4f).

The first example is a GPR profile collected over an aeolian dry dune in the Chadian desert (Bano et al., 1999). The goal of this survey was to image the interior and the base of the dunes to better understand sedimentological processes. The GPR profile has been obtained using a 450-MHz shielded antenna. The acquisition mode was a constant offset of 0.25 m, the antennas were moved by 0.125 m steps with a stack of 64 to improve the signal-to-noise ratio.

A standard processing (with in-house interactive GPR software) has been applied and the resulting profile is displayed in Fig. 5a. The following processing sequence was used: constant shift to adjust the time zero followed by normal move-out corrections; running average (DC) filter to remove the low frequency; flat reflections filter to remove some clutter noise (continuous flat reflections) caused by multiple reflections between shielded antennas and the ground surface; a band-pass filter and finally a time-varying gain function. The same standard processing is applied to all GPR data presented in this section.

The GPR profile of Fig. 5a shows complex geometries, with imbricate reflections corresponding to different deposit phases. The undulating reflection indicated by four white arrows in this figure represents the base of the dune, which in fact is nearly flat and consists of pebbles (> 2.0 mm in diameter). This reflection stems from the contact between the aeolian sands near the surface and deeper lake deposits consisting of an unconsolidated silty sandstone layer of very fine to medium grain size. In order to apply the topographic static shift and/or migration, we need to know the velocity of the GPR waves. In Fig. 5a, we also observe a nice 10-m-wide (80 traces) diffraction hyperbola situated just under the base of the dune (see black circle). After analyzing this diffraction, with different velocities, we found that it can be fitted very well with a constant velocity model of 0.18 m/ns. This value represents an average velocity from the surface of the dune to the diffraction point and it is in good agreement with values found in the literature for dry sands (Gómez-Ortiz et al., 2009; Guillemoteau et al., 2012).

Fig. 5b shows the same GPR profile as in Fig. 5a, but when a standard migration followed by topographic corrections is performed, using a velocity of 0.18 m/ns. The topography shows a local variation of about 30% (at profile coordinate 38 m, black arrow in Fig. 5b) and its global variation of about 5 m (5%) is comparable to the investigation depth. The diffraction hyperbola is well collapsed (at profile coordinate 25 m, black circle) and the reflection from the base of the dune is roughly flattened. Below the area of high topographic gradient (38 m horizontally), we observe a very bad feature (black arrow). The whole area looks blurred, and reflectors are losing consistency. The results of the standard migration followed by topographic corrections are bad.

Fig. 5c presents the profile after a static shift followed by a standard migration. The migration hyperbola (black circle) is slightly over-migrated. The bad feature indicated by the black arrow in Fig. 5b is corrected. The reflectors are now consistent and the dipping reflector shown by the black arrow (Fig. 5c) has been moved up-dip.

Fig. 5d presents the topographic migration with the same velocity (0.18 m/ns, as in both previous cases) and a specific migration template 13 m wide (100 traces) has been chosen, which is slightly larger than the width of the observed hyperbola (10 m) on the profile. The base of the dune is flattened and the diffraction at 25 m is now correctly focused on a single point inside the black circle, which justifies our choice of 0.18 m/ns for the GPR velocity. The dipping reflector shown by the black arrow has undergone a vertical and horizontal shift of 1.1 and 3.8 m, respectively. It starts at the base of the dune and goes up-dip rightwards as expected (on the non-migrated section of Fig. 5a, these reflections were crossing the base of the dune).

The measured dips on the topographic migrated section of the same reflectors (shown by white arrows) are slightly larger than the dips measured in Fig. 5c (static shift followed by standard migration). Their values are now 26.5° and 19.5° on the topographic migrated section, instead of 25° and 17.7° in Fig. 5c. Although the global topographic variation of the profile does not exceed 5%, the result of the topographic migration is slightly better than the result of the static shift followed by the standard migration. Remember here that the later routine over-migrates the data at large depth (case of the diffraction under the base of the dune).

Fig. 6 shows a portion of the profile of Fig. 5d with topographic migration for different velocities ranging from 0.16 m/ns (left) to 0.20 m/ns (right) with a 0.1 m/ns increment. The hyperbola is not collapsed for the two first figures, while it is over-migrated for the last two ones. The middle figure shows the migrated image with the correct velocity of 0.18 m/ns. The depths and the dips of the reflectors are also changed. The depth of the diffracting point is ranging from 6.5 m (for a velocity of 0.16 m/ns) to 7.8 m (for a velocity of 0.20 m/ns). Therefore, a change of around 5% in the velocity causes a change in depth of nearly 0.3 m (for a depth of around 7 m). The dip of the reflector indicated by the arrow is ranging from 22.6° (velocity of 0.16 m/ns) to 29.9° (velocity of 0.20 m/ns). The dip increases by roughly 2° per 0.01 m/ns velocity increase. To have a correct migration, we conclude that the precision in the estimation of the velocity should be better than 5% of the true velocity.

In 2010, we conducted a GPR campaign in Mongolia, 80 km to the west of the capital Ulaanbaatar. The context of this study was seismic hazard. Fig. 7 shows a GPR profile obtained with an unshielded 50 MHz Rough Terrain Antenna (RTA). The profile is more than 200 m long and is perpendicular to the Hustay Fault. This is in a context of a very low slip rate (most likely less than 1 mm per year), and the fault geomorphology has been smoothed during a long period of erosion. Therefore displacements in the topography are not observable. However, in the field, evidences of the fault plane are still visible. Most of the profiles acquired in this area display a strong reflection, which corresponds to the fault plane. These profiles give complementary information such as the dip of the structure and the exact location of the fault near the surface to help design the layout of future palaeoseismic campaigns. The acquisition mode was a constant offset of 4.1 m; traces have been recorded every 0.2 m with a stack of 16 to improve signal-to-noise ratio.

The processing used to obtain Fig. 7a is similar to the one used in the case of the Chadian GPR data. A velocity analysis, which is not presented here, has been done over the surveying area by analysing diffraction hyperbolae present in the GPR data. A mean velocity of 0.12 m/ns has been determined for the whole area. As in the previous example, the topographic variation of 20 m (slope less than 10%) is comparable to the investigation depth. Figs. 7b, c and d respectively display the data after standard migration followed by static shift, static shift followed by standard migration, and topographic migration. The diffraction hyperbola indicated by the black circle is well focused in Fig. 7b and d, and appears slightly over-migrated in Fig. 7c (as in the case of the Chad dune). The dipping reflector (fault plane) indicated by the white arrow now displays a constant slope down to a depth of 24 m in Fig. 7b and d. However, on the section of Fig. 7c, the reflector is attenuated at a depth of 17 m and has lost its continuity. Its dip angle is changing from 32.2° (migration and static shift) to 34° (static shift and migration and migration with topography). The main observation is the location of the reflector, which is very similar in the cases of Fig. 7b and d, while in Fig. 7c the reflector has been shifted (5.5 m horizontally and 2.6 m vertically) and is reaching the surface. In this case, the migration followed by static shift seems to give more convincing results than the static shift followed by migration (which is the opposite of what was observed in the Chadian dune example). From this result, we conclude that topographic migration should be considered at any location where the subsurface shows steep dip angle structures (exceeding 30°), even in the case where the surface slope is less than 10%.

As in the previous case, we performed a velocity sensitivity analysis by using a topographic migration of the fault plane reflection with different velocities ranging from 0.1 m/ns to 0.14 m/ns with a 0.01 m/ns step. After topographic migrations, the slope of this reflector is varying from 28.9° to 39.5° and increases by roughly 2.65° for a 0.01-m/ns increase in the velocity. In this case, for a correct topographic migration, the estimation of the velocity should be better than 8%.