New geophysical technologies have been developed to solve the specific geological survey problems in areas where it is hard to access on the ground, such as mountainous areas with acute topographic relief, lakes and bogs(Yang et al., 2013; Ge et al., 2009). Of them, the ground-airborne transient electromagnetic method is a useful tool. It places galvanic or loop sources on the ground to excite the high-power transient electromagnetic field and uses unmanned aerial vehicle or airship with a probe to collect data in the air. Compared with the airborne transient electromagnetic method, the data of this new method has higher signal-to-noise ratio due to ground-based sources used. In addition, measuring with an unmanned aerial vehicle or airship simply avoids certain safety issues of a pilot. Meanwhile, this method can also highly improve work efficiency, especially in those aforementioned areas.

The development and application of the ground-airborne transient electromagnetic method will bring revolutionary change to the traditional ground transient electromagnetic methods. Denser survey points can help realize three-dimensional survey coverage. It can adopt single sources or multiple sources. And the pseudoseismic imaging method can be used for its interpretation. With more sources and the higher-power exciting field, investigation depth of the ground-airborne transient electromagnetic method can be increased by several kilometers. And on this basis, the transient electromagnetic method can be introduced to investigation toward depth and used for three- dimensional geological mapping.

The idea of the ground-airborne transient electromagnetic method was first proposed by Nabighian in 1988, in which an electric dipole as its transmitting source is used. Mogi et al.(1998, 2009)developed the GREATEM system and managed to apply it to studying the characteristics of sea/l and boundaries in coastal areas and also in a volcanic structure survey. Smith et al.(2001)compared the survey data of airborne, semi-airborne and ground transient electromagnetic methods. Verma et al.(2010)analyzed the response characteristics of the GREATEM system on a half-space model. Yang et al.(2012)studied the data preprocessing of ground-airborne electromagnetic methods in the galvanic source time domain. Li et al.(2013)studied the ground-airborne electromagnetic signals de-noising using a combined wavelet transform algorithm. Until now, a complete ground-airborne transient electromagnetic system has not been established due to the lack of three-dimensional interpretation methods.

The research issues relating to ground-airborne transient electromagnetic interpretation include the following parts: full field apparent resistivity definition, Kirchhoff migration of virtual wave field, and reverse synthetic aperture imaging of the ground-airborne transient electromagnetic method. Based on existing research achievements, this paper builds a complete ground-airborne transient electromagnetic interpretation system with three-dimensional TEM pseudo-seismic imaging and lays the foundation for the entire investigation system.

The full field apparent resistivity definition is a key problem in ground-airborne electromagnetic investigation. In practice, it is desirable to acquire apparent resistivity without limitation in distances, time periods and flight heights. With such full-field apparent resistivity, a ground-airborne investigation can be implemented in three-dimensional data collection. With three-dimensional data processing and three-dimensional interpretation, observation technologies and interpretation of the transient electromagnetic method can be improved to a higher level.

As shown in Fig. 1, source AB is a horizontal electric dipole placed on the surface of the homogeneous earth in X direction, with a length of ds and current I = I₀e^-iωt. The origin of coordinates O is the central point of AB. The Z axis points downwards. The projection on the surface of any point M in the air is P. The distance from O to P is r, and the angle between r and X direction is φ.

(Fig.1)

Fig.1 Diagram of an electric dipole on the surface of the homogeneous earth

From the Helmholtz equation, with boundary conditions on the earth surface of electromagnetic field, the curl relation between magnetic field and electric vector potential, we can get the vertical component of magnetic field in the frequency domain at any point in the air produced by an electric dipole on the homogeneous earth:

(Fig.2)

Fig.2 Vertical component of the magnetic field at different times with different offsets in homogeneous half-space (a)Offset 50 m;(b)Offset 500 m;(c)Offset 1500 m;(d)Offset 3000 m.

Giving an initial value and Taylor exp and H_z(ρ, r, t)in the neighborhood of ρ₀:

Figure 3 is the full field apparent resistivity of two-layered models with various resistivity values of the second layer, different offsets and several flight heights. The length of the galvanic source is 1000 m and the rest of the model parameters are ρ₁ = 100 Ωm, h₁=20 m, ρ₂=2, 5, 10, 30, 80 Ωm. No matter at high flight height of 100 m or low flight height of 50 m, a long offset of 3000 m or short offset of 500 m, the full field apparent resistivity can always converge to the values of first and last layers in the model at early and late times, respectively. It can also show the electrical changes at transitional times.

(Fig.3)

Fig.3 Full field apparent resistivity of two-layered models with various resistivities of the second layer at both short and long offsets at different receiving heights (a)Receiving height 50 m, offset 500 m;(b)Receiving height 50 m, offset 3000 m;(c)Receiving height 100 m, offset 500 m;(d)Receiving height 100 m, offset 3000 m.

Figure 4 is the full field apparent resistivity of three-layered models with various resistivity values of the second layer. Model parameters are: H model: ρ₁=100 Ωm, ρ₂=5, 10, 30, 50 Ωm, ρ₃=100 Ωm, h₁=20 m, h₂=10m; K model: ρ₁=10 Ωm, ρ₂=30, 60, 100, 200 Ωm, ρ₃=10 Ωm, h₁=10 m, h₂=20 m; A model: ρ₁=30 Ωm, ρ₂=50, 60, 70, 80 Ωm, ρ₃=100 Ωm, h₁=10 m, h₂=20 m; Q model: ρ₁=100 Ωm, ρ₂=20, 30, 50, 80 Ωm, ρ₃=10 Ωm, h₁=20 m, h₂=10 m. The first and last parts of all the full field resistivity curves can reach ρ₁ and ρ₃, respectively. The curves have obvious differentiation with parameter change. Compared with apparent resistivity derived from ∂B(t)/∂t, those derived from with B(t)can better identify electrical change.

(Fig.4)

Fig.4 Full field apparent resistivity of three-layered models with different resistivities of the second layer at two typical offsets (a)H model, offset 500 m;(b)K model, offset 500 m;(c)A model, offset 500 m;(d)Q model, offset 500 m;(e)H model, offset 3000 m;(f)K model, offset 3000 m;(g)A model, offset 3000 m;(h)Q model, offset 3000 m.

As Kirchhoff integration immigration uses the wave field during processing, we should first transform the transient electromagnetic field into a virtual wave field(Lee et al, 1989; Hoop et al, 1996; Chen et al, 1999; Li et al, 2001, 2005, 2011, 2013). The transient electromagnetic field and virtual wave field have the following relationship:

With the virtual wave field and continuous velocity distribution of underground space, we can get the Kirchhoff integration solution from the wave equation(Li et al., 2010):

Immigration is the inverse process of data recording. Setting upward waves of the wave field as G $\left( \xi ,\eta ,{{\varsigma }_{0}},t+\frac{r}{v} \right)$, Eq.(10)then can be rewritten as

Applying the three-dimensional boundary element method to Eq.(12)can implement the immigration of the virtual wave field. This is actually the inverse process of data recording. Then the locations of virtual wave sources of underground reflecting interfaces and the spatial distribution of underground geological bodies can be determined.

Synthetic aperture imaging in airborne transient electromagnetic methods are resemble to the synthetic aperture radar system(Li et al, 2012; Qi, 2012). It composites smaller sized real antenna apertures to obtain a larger sized transmitter loop with equivalent aperture by means of data processing to improve the resolution and the detection depth.

An inverse synthetic aperture radar system usually uses a fixed radar to get the longitudinal and horizontal imaging with a high resolution of a moving object, which seems opposite to the synthetic aperture radar system but is actually the same in essential(Berrizzi et al, 1996; Bao, et al, 2001; Li, 2011). The transmitter of the ground-airborne transient electromagnetic method is fixed on the ground while the receiving device is in the air, which is similar to the inverse synthetic aperture radar system. With the theory of the inverse synthetic aperture radar system, we apply correlation stack processing to signals with a certain range around the survey points in the air to get a larger sized receiving loop to further improve the longitudinal and horizontal resolutions of detection.

Figure 5 is a sketch of inverse synthetic radar imaging, where A and B indicate the two current electrodes on the ground which excite the transient electromagnetic field in the subsurface. When the transient electromagnetic field reaches the geological body, it will induce secondary eddy-currents. After the transmitter current is cut off, electromagnetic signals decaying with tine of the secondary magnetic field induced by a secondary eddy-current can be received in the air. As the electromagnetic signals within the L aperture range from the same geological body have correlation characteristics, we take the central point of the L aperture range as the reference point i and mark the wave field of this point as U(r_i, τ), here r_i is the distance from point i to any point(with a serial number of -N, · · ·, N)within the L aperture range, and τ is the relative moving amount. We take the maximum correlation coefficient as the weight coefficient, which is derived from correlation calculation of the reference point with each point within the L aperture range. We multiply the weight coefficient with the wave field value of each point and stack them to the central point, the composite value of the central point can be written as

(Fig.5)

Fig.5 Sketch of inverse synthetic aperture imaging

A model with a mined out water goaf is designed to test the efficiency of our method. As shown in Fig. 6, the goaf is placed within a coal bed. The upper and middle layers of the coal bed are fine s and stone and arenaceous s and stone, respectively. The lower layer of the model is quartz s and stone. Parameters of the model are listed in Table 1.

(Fig.6)

Fig.6 Diagram of a complex model with one mined out water body

(Table 1)

Table 1 Parameters of a complex model with one mined out water body

Table 1 Parameters of a complex model with one mined out water body

Figure 7 is the top view of the model. With the three-dimensional diagram in Fig. 6, we can get the location of the coal bed and the goaf. In total 301×301×100 cubic uniform grids with a side length of 10 m are used in forward modeling. The length of the galvanic source AB is 100 m. The central point of AB is set at the center of the model. Three survey lines of LineX81, LineX86 and LineX91 with a flight height of 100 m are picked out for data processing, whose offsets are 0 m, 50 m and 100 m, respectively(see Fig. 7). The full field apparent resistivity contour maps of these three selected survey lines are plotted in Fig. 8. As LineX81 crosses the central region of the goaf, the low-value apparent resistivity contours clearly indicate the location of the goaf. LineX86 is around the boundary of the goaf, but the apparent resistivity contours can still show the existence of the underground low resistivity body. While LineX91 is far away from the goaf, the anomalies of the goaf disappear. In total, the full field apparent resistivity contours derived from our method coincide well the distribution of the goaf.

(Fig.7)

Fig.7 Top view of the complex model with one mined out water bodies and position of the transmitting source and the survey line arrangement

(Fig.8)

Fig.8 Apparent resistivity contours of the selected lines (a)Apparent resistivity contour of LineX81;(b)Apparent resistivtiy contour of LineX86;(c)Apparent resistivity contour of LineX91.

We apply wave field transform to the electromagnetic field with a flight height of 100 m from the threedimensional forward modeling(Fig. 9). The results show that there are two interfaces, with the upper one indicating the ground-air interface. As the conductivity contrast between the air and the earth is large, the wave field signals are high in amplitude and spread all over the region. While the strong signals of the lower interface focus in the central part and become weak in surrounding areas, it obviously indicates the existence of the goaf.

(Fig.9)

Fig.9 Results of 3-D wave field transform:(a)Three-dimensional diagram;(b)Slice diagram

The reverse synthetic aperture imaging technique can be applied before or after Kirchhoff immigration. Fig. 10 is the result of an immigration from the wave field transform and Fig. 11 is the result of correlation from migration. It can be seen that the ground interface changes slightly after correlation, while the goaf region narrows with clear boundaries. The effect of reverse synthetic aperture imaging is good.

(Fig.10)

Fig.10 Results of 3-D migration imaging: (a)Three dimensional diagram;(b)Slice diagram

(Fig.11)

Fig.11 Results of inverse synthetic aperture imaging (a) Three-dimensional diagram; (b) Slice diagram; (c) Diagram after removing high-energy reflection.

(1)Defining full field apparent resistivity based on reverse function and iteration algorithm makes it possible to calculate apparent resistivity in full time and full space. The calculation results of layered earth models using this method show that the curves of resistivity reveal obvious differentiation of earth electrical change.

(2)The synthetic models show that apparent resistivity from the electromagnetic field has no false anomalies. The first and last part of the apparent resistivity curves can approximate the uppermost and lowest layer of the model, respectively.

(3)The pseudo seismic algorithm for the transient electromagnetic method can help realize three-dimensional imaging of the transient electromagnetic field and improve the interpretation to a higher level.

(4)The reverse synthetic aperture imaging for the ground-airborne transient electromagnetic method is proposed to improve the resolution with a reverse synthetic aperture radar system. The synthetic model with a goaf proved that correlation stack can strengthen useful signals, improve the signal-to-noise ratio and the resolution, which verifies the effect of reverse synthetic aperture imaging for the ground-airborne transient electromagnetic method.

This work was supported by the National Natural Science Foundation of China(41174108).