Passive Super-Low Frequency electromagnetic prospecting technique

Abstract

The Super-Low Frequency (SLF) electromagnetic prospecting technique, adopted as a non-imaging remote sensing tool for depth sounding, is systematically proposed for subsurface geological survey. In this paper, we propose and theoretically illustrate natural source magnetic amplitudes as SLF responses for the first step. In order to directly calculate multi-dimensional theoretical SLF responses, modeling algorithms were developed and evaluated using the finite difference method. The theoretical results of three-dimensional (3-D) models show that the average normalized SLF magnetic amplitude responses were numerically stable and appropriate for practical interpretation. To explore the depth resolution, three-layer models were configured. The modeling results prove that the SLF technique is more sensitive to conductive objective layers than high resistive ones, with the SLF responses of conductive objective layers obviously showing uprising amplitudes in the low frequency range. Afterwards, we proposed an improved Frequency-Depth transformation based on Bostick inversion to realize the depth sounding by empirically adjusting two parameters. The SLF technique has already been successfully applied in geothermal exploration and coalbed methane (CBM) reservoir interpretation, which demonstrates that the proposed methodology is effective in revealing low resistive distributions. Furthermore, it siginificantly contributes to reservoir identification with electromagnetic radiation anomaly extraction. Meanwhile, the SLF interpretation results are in accordance with dynamic production status of CBM reservoirs, which means it could provide an economical, convenient and promising method for exploring and monitoring subsurface geo-objects.

Keywords

Super-Low Frequency (SLF) / three-dimensional modeling / frequency-depth transformation / geothermal exploration / coalbed methane (CBM) / electromagnetic radiation (EMR)

Introduction

Conventional remote sensing techniques, with respect to regular bands (visible/infrared/microwave bands), are commonly used in land surface imaging (Atzberger, 2013; Congalton et al., 2014) and near-surface surveys (Yavuz and Teixeira, 2009; Slonecker et al., 2010). However, these techniques utilize electromagnetic (EM) waves which are useless in the subsurface exploration of tens of meters due to the strong attenuation in the near-surface earth (Ouchi, 2013). Advanced remotely sensed systems have recently been developed to increase the observation ranges (Melesse et al., 2007; van Genderen, 2011). As an example, prospecting devices were created that use natural source EM signals at a wide frequency, ranging from 3–3000 Hz (Wang et al., 2014b), covering the Extremely-Low Frequency (ELF, 3–30 Hz), the Super-Low Frequency (SLF, 30–300 Hz), and the Ultra-Low Frequency (ULF, 300–3000 Hz). On the basis of model tests, the SLF range is the most significant frequency range applied in the subsurface exploration, and is explicitly defined by ITU. Thus, the depth sounding technique using the low frequency range (3–3000 Hz) is named after the generalized “SLF.” Furthermore, as a complimentary remote sensing tool, the SLF electromagnetic prospecting technique (SLF technique) has undergone a revolution which is driven by its increasing use in subsurface observation. It has become one of the most applicable remote sensing methods for subsurface imaging, especially for geological exploration (Li et al., 2013; Wang et al., 2014c). According to the “skin depth” formula in the classical EM theory, SLF EM waves can penetrate a few kilometers below the surface while the frequency decreases. The corresponding responses provide substantial geological information which can be observed at the surface (Qin et al., 2010). Therefore, the SLF method increases the depth range of the earth observation by as much as 10 km showing a significant advantage over regular wavebands, including lower attenuation, solid theoretical support, and all-weather adaptivity. These wavebands are more in favor of underground geo-object detection. During the past decade, hardware and software systems of the SLF technique have been designed and improved (Cui et al., 2009; Shan et al., 2014). The data processing algorithms including independent component analysis, the wavelet and lifting-wavelet analysis are being adopted (Wang et al., 2013). Some preliminary interpretation techniques are now applied in many geological aspects (Qin et al., 2010; Wang et al., 2014b).The feasibility of the SLF technique has been theoretically and practically explored in coalbed methane (CBM) reservoir interpretation on the basis of numerical modeling of 2-D reservoir models (Wang et al., 2014b). Nevertheless, previous studies were not able to provide a complete illustration of theoretical basis, depth resolution, and frequency-depth transformation, nor provide a fully systematic application scheme. Due to the 3-D characteristics of Earth, theoretical analysis of realistic 3-D models for the SLF technique still needs to be fully explored. When taking the convenience and accuracy of finite difference methods into consideration, it is clear that 3-D numerical algorithms, based on the staggered finite difference grids, are becoming one of the most adopted methods for calculating natural source EM responses of realistic 3-D geo-electrical models (Mackie et al., 1993; Siripunvaraporn et al., 2002). In addition, regular computers or workstations have the capability to accommodate the computational costs of 3-D modeling (Börner, 2010).

In this paper, the basic theory and numerical algorithms of the SLF technique are discussed. Theoretical SLF responses of multi-dimensional models are then calculated, and feasibility analysis, including the depth resolution and frequency-depth transformation are further discussed in detail. Consequently, a complete work flow and case studies are explored. This study demonstrates that the SLF technique would be a beneficial subsurface geological sensing tool in the near future.

Super-Low Frequency electromagnetic prospecting theory

Basic theory

The SLF EM prospecting method originated from the magnetotellurics (MT) method. Both methods use the natural source EM signals induced within the Earth and measurements of EM field data with transformations into the frequency domain, which allow us to map the subsurface geo-electrical distribution as a function of position and depth (Meju, 2002). However, there is a difference in that the MT method uses the impedance tensor and tipper as response functions which are estimated from EM data in multiple channels (i.e., electrical and magnetic channels). In contrast, we propose a single horizontal magnetic amplitude component as the SLF response (Wang et al., 2014b). Since the horizontal magnetic amplitude from a single site contains resistivity information about the Earth’s structure, the proposed methodology is straightforward in the implementation of the one-dimensional (1-D) models. Starting with the governing Maxwell equations in the frequency domain:

(1)\nabla \times \overrightarrow{H}=\overrightarrow{E}/\rho ,

(2)\nabla \times \overrightarrow{E}=i\omega \mu \overrightarrow{H},

where \rho is the resistivity, \omega is the angular frequency/Hz, and \mu is defined as the magnetic permeability. Here, \mu is approximated as the permeability in the air according to the MT theory, and only the resistivity is a spatial variable in the 3-D conducting medium. Eqs. (1) and (2) can yield the Helmholtz equations based on the time-harmonic electrical field:

(3)\nabla \times \nabla \times \overrightarrow{E}=i\omega \mu \overrightarrow{E}/\rho (x,y,z),

or they can give the Helmholtz equations based on the time-harmonic magnetic field:

(4)\nabla \times \rho (x,y,z)\nabla \times \overrightarrow{H}=i\omega \mu \overrightarrow{H},

where \rho (x,y,z)in Eqs. (3) and (4) represents the complex-valued resistivity at different spatial positions.

A Cartesian coordinate system can be defined with z downward as the vertical direction, and assumes that \rho (x,y,z) varies only with z direction. The resistivity \rho (x,y,z) is invariant in the horizontal x and y directions. Such simple models are deemed as 1-D models, and especially for the half-space, \rho is the homogeneous resistivity. Eq. (3) can be simplified as:

(5){\displaystyle \frac{\partial E_x^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}-k^{\mathrm{2}}{E}_x=\mathrm{0},

(6){\displaystyle \frac{\partial E_x}{\partial z}}=i\omega \mu H_y,

where k=\sqrt{-i\omega \mu /\rho }. Horizontal magnetic amplitude can then be obtained from Eq. (6):

(7)H_y(z)={\displaystyle \frac{\mathrm{1}}{\sqrt{-i\omega \mu \rho }}}A{e}^{-kz},

where A is a constant. The SLF magnetic amplitude correlates with the resistivity \rho distribution of geo-objects. Consider a layered model with N stacked layers beneath the surface. It is easy to derive a recursive expression of the SLF magnetic component:

(8)H_y(z)=-{\displaystyle \frac{k_m}{i\omega u}}{A}_m\left(e^{-k_{m}z}-{\displaystyle \frac{B_m}{A_m}}{e}^{k_{m}z}\right),

where k_m=\sqrt{-i\omega \mu /{\rho }_m}, {\rho }_m is the resistivity of m^th layer. A_mand B_m are both constants, and more specifically, {\displaystyle \frac{B_m}{A_m}} can be solved by the iterative impedance solution. Letting z=0 (i.e., the acquisition site is located on Earth’s surface), the normalized amplitude of H_y(z) is explicitly a stable SLF response which only depends on the geo-electrical structure beneath the surface.

According to the MT theory, any linearized source can be divided into the combination of two orthogonal sources in XY and YX modes (Avdeev, 2005). Earth’s multi-dimensional conducting subsurface is revealed through assuming itself as a 1-D layered model, 2-D models in the TE and TM modes (i.e., XY and YX modes), as well as 3-D models in XY and YX modes. With the aim of capturing the SLF response characteristics of most geologically plausible structures, numerical methods are required to solve the forward calculations of the SLF magnetic amplitude responses. In this regard, only 1-D models have analytical solutions of normalized SLF responses. 2-D and 3-D models require sophisticated numerical finite-difference, finite-element, or integral equation methods to estimate the SLF response characteristics. The developed algorithms are all based on different equations and discretized meshes, which can accurately determine 3-D EM response values of different nodes on the surface (Siripunvaraporn et al., 2002). As the most common manifestation of discretization methods, the finite-difference method is simple to construct and economically computational, and has hence been adopted in this paper. The 3-D model is fully solved using an approach facilitated by staggered grids (Mackie et al., 1993). The key step for setting up the large-scale linear equations is to first implement the boundary conditions on the perimeter of the computational area. Setting the electrical field value equal to a constant on the upper boundary is reasonable. The bottom boundary should approximate the 1-D solution for consistency with the stitched 1-D assumption. Subsequently, the 2-D numerical solutions are set on the left and right sides (Aprea et al., 1997). The entire computational domain can be discretized using the integral forms of Eqs. (3) or (4). The solutions produced from the electrical field values in Eq. (3) are closer to the exact solutions with a high numerical accuracy and easy convergence (Siripunvaraporn et al., 2002). As a result, our algorithm has been made suitable for the calculation of SLF responses. The large-scale, sparse and complex linear equations, based on electrical field values, are then solved using a stable quasi-minimal-residual iterative solver with a Jacobi pre-conditioner. The interpolation method is also used to calculate the approximate value of the magnetic field using Eq. (2). Consequently, the theoretical magnetic amplitude spectrum are also obtained (Börner, 2010). Finally, the SLF responses and MT responses in the frequency domain are calculated.

In order to undertake an analysis of computational accuracy of the algorithms, a simple 1-D model was designed. Figure 1(a) shows a 1-D layered model consisting of 1W·m 30 m thick layer surrounded by 100 W·m layers over a 3000 W·m basement. The basic aim of this model yields analytical MT apparent resistivity responses for examining the performances of 1-D, 2-D, and 3-D algorithms. In view of the 1-D structure, horizontal grid discretization with regard to 2D and 3D algorithms should have no effect on MT responses. Thus, the acquisition site is located at the model surface center. Figure 1(b) shows few discrepancies of apparent resistivity responses at small scale for 2D & 3D numerical algorithms in comparison to the 1-D analytical solutions.

Furthermore, algorithms were verified using general 3-D geo-electrical models proposed by the Comparison Of Modeling Methods for ElectroMagnetic Induction (COMMEMI) and adopted by academic communities (Zhdanov et al., 1997). One representative model configuration of COMMEMI is Model 3D-1. This model is a rectangular 0.5 W·m conductive anomaly with dimensions of 1 km×2 km×2 km when inserted in a homogenous 100 W·m half-space. Calculations for this model were carried out along two coordinate axes for the anomaly detection at the Earth’s surface for one frequency point (10 Hz). Even though the 3D-1 model is geometrically simple, a high contrast of conductivity can result in a relatively large anomalous EM field near the surface. The mesh discretization for the 2-D algorithm requires 24 blocks in the horizontal direction and 29 blocks in the vertical direction (with 7 air layers). Alternatively, the length along another horizontal direction for the 3-D algorithm was finely divided into 22 blocks. The other discretization was the same with the 2-D configuration. A complete comparison of the results are apparent resistivitys in XY and YX modes calculated by Russian, American, and German scientists using integral equations as shown in Fig. 1(c). Our proposed 2-D and 3-D algorithms have been shown to be computationaly accurate in comparison to the recommended results by COMMEMI in both XY and YX modes. This demonstrates that 2-D and 3-D algorithms are reliable for further use in the calculation of SLF responses.

Theoretical SLF responses

The SLF theoretical responses of the model in Fig. 1(a) using aforementioned 1-D, 2-D, and 3-D numerical algorithms are depicted in Fig. 2. The results of the horizontal magnetic components in some modes are similar and thus are not given. By taking into consideration the variance in absolute values of SLF magnetic responses for such factors as boundary conditions, meshes, and differential approximations, all the solutions were normalized between 0.1 and 1. The resultant theoretical responses are nearly equivalent, and reveal that the SLF responses show a tendency toward power laws. Furthermore, the accuracy of our proposed algorithms are acceptable for simulating SLF responses.

In geological survey, near surface electrical anomalies at local scales can be deemed as 3-D anomalies compared to the regional 2-D distributions of electrical conductivity. The objective layers are saliently generalized as 2-D geological structures. To illustrate 3-D distortions in SLF responses, a 3-D model, combined with the high resistive basement layers and relatively uniform adjacent layers, is constructed and simulated using developed finite difference algorithms. This model is designed matching the exploration ranges of the SLF method in Fig. 3(a): in the XOZ profile (Y=0), Y is the strike direction. The model consists of a layered background impregnated with a conductive body at the depth of 500 m. The conductive target is 3000 m× 30 m with 1 W·m. The shallow part contains 3-D resistive (1000 W·m) and conductive (1 W·m) anomalies (1000 m on the x axis, 50 m on the z axis). The background spanning from surface to 600 m in depth is 100 W·m, over a 3000 W·m semi-infinite basement layer. The XOY profile (Z=0) is depicted in Fig. 3(b), and the length of two anomalies on the y axis is 1000 m. The objective layers extend infinitely along the strike direction to a 2-D structure. Two typical sites (shown by black points) were selected to analyze SLF magnetic amplitude responses. The model was discretized into 32×24×28 grids for 3-D modeling. The anomalous bodies should be carefully subdivided to obtain more accurately simulated SLF magnetic responses.

Figures 3(c) and 3(d) illustrate the normalized SLF magnetic responses and correspondingly average responses of two sites. Notice that the shallow 3-D anomalies dominate the discrepancy in the magnetic amplitude responses in different modes or directions. Previous research demonstrates that near surface electrical anomalies could lead to a static shift owing to which galvanic amplification or reduction in electrical field signals is superimposed on 2-D electrical responses (Berdichevsky et al., 2010). The 2-D responses are only obtained from 2-D models without near surface anomalies. This contrast indirectly reveals the 3-D distortions on the 2-D magnetic amplitudes, and the impacts when large or small variations are placed on magnetic amplitude signal levels in different modes and directions. If measured, the average response of all modes and directions would more accurately match the 2-D responses. This same result is observed forthe site above high resistive near-surface anomalies as shown in Fig. 3 (d). As for the SLF responses of the model with a near-surface conductive anomaly, the average response still shows some deviation even though it is closer to the exact 2-D response than the 3-D response. To reduce the 3-D distortions in SLF responses, an average response using more repitition of acquired data for one site has been performed through model tests. This indicates that the reasonable SLF responses for interpreting objective geo-targets should use average normalized responses over repetitive data (i.e., total 10 sets of data in different modes and directions for each site).

Feasibility analysis

Feasibility analysis is the key factor for determining the potential of the SLF technique. On the basis of the MT theory, the SLF magnetic amplitude changes correlatively with the resistivity distributions of geo-objects. For most plausible sedimentary structures, hypothetical resistivity models, rather than simplicity models, are intended to be an idealized representation of three-layer models, i.e., cover, objective layer, and basement. It follows then that feasibility analysis is based on the response analysis of these hypothetical models (i.e., “objective layer” models), including evaluations of the depth resolution and Frequency-Depth transformation. Resistivity contrasts and objective layer scales should be implemented, which are subjective in the following theoretical analysis.

Depth resolution of objective layers

In order to explore the depth resolution of objective layers, the impacts of objective resistivity and the ratios of “objective depth to layer thickness” on SLF responses were numerically evaluated on the basis of the previous research of Huang and Won (2003). With the aim of conducting geological surveys at a depth of 1000, we designed a series of three-layer models, including the cap layer, the objective layer, and the basement layer. The cap layer is 100 W·m with various thicknesses, and the objective layer is 100 m thick with different resistivity values over a 1000 W·m uniform basement layer in a semi-infinite half space.

Figure 4(a) depicts the responses of models when the cap layer is 500 m thick and the resistivity values of objective layers change from 0.01 W·m to 10000 W·m. When the objective layer is more conductive, especially by two orders of magnitude, the response curves in the frequency domain rise at low frequency compared to the gently changing trend line (dotted line), which is manually drawn as an approximate straight line. The uprising amplitudes are more obvious while the objective layer is more conductive. On the other hand, when the frequency is extremely low, the SLF EM waves penetrate as deep as the high resistive uniform basement layer. As a result, the SLF responses are stabilized. If the middle layer is 10 W·m, one order of amplitude is more conductive than the adjacent layers, then the uprising tendency of the SLF responses at low frequency is nearly undetectable. However, the downward tail at extremely low frequency still approaches a stable value. When the resistivity values of middle layers exceed 10 W·m, the anomalous responses from objective layers cannot be identified. Thus, the SLF responses have proven to be suitable for resolving conductive anomalies. If the deep objective layers are one or more orders of amplitude, i.e., more conductive than the adjacent layers, then they can be clearly interpreted. However, resistive layers cannot be distinguished from the background.

When the objective layer is set to 10,000 W·m, and the depth ranges from 200 m to 1000 m, the “depth to thickness” ratios range from 2 to 10. A series of SLF responses are depicted in Fig. 4(b). When the middle layer is deeper than 800 m, the frequency point (i.e., the turning point), where the curve turned downward in contrast with the trend line, is much lower with a slower variance. Therefore, the resultant SLF responses of these models tend to be the same. Thus, when the middle layer is 100 m and the ratio of “the objective depth to its thickness” is over 8, the turning points are almost the same.,As a result, the depths of the resistive objective layers could not be identified. A series of models for the middle conductive layers of 0.1 W·m, with objective depths from 200 m to 1000 m, were also established. The results are shown in Fig. 4(c). The turning frequency points are lower if the middle layers are deeper. At lower frequency, a stable value is approximated in the tails of SLF curves. When the ratios of “depth to thickness” are over 8, the objective depths are difficult to interpret, even though the SLF response features of conductive objective layers are clearly compared to dotted trend lines. Therefore, the SLF curves can be reasonably used in the direct interpretation of geological structures.

Frequency-depth transformation

Another available key component of the SLF technique allows for the immediate estimation of the underlying resistivity profile through a parametric representation of a frequency-depth transformation. The SLF responses in the frequency domain are transformed into an approximate depth profile with the frequency as an independent parameter. This kind of study is mainly available as the direct inversion in the MT theory. The first considered method is the \rho *(z*) method (Schmucker, 1987). The physical basis of the method is based on the identification of Schmucker’s response as the normalized first moment of induced in-phase currents. However, it is ineffective in common geological instances with a thick resistive layer below a conductor. The other classical direction inversion method is Bostick inversion (Gomez-Treviño, 1996; Rodríguez et al., 2010). Bostick inversion attains the asymptotic lines of apparent resistivity response curves in a double-logarithmic plot when the basement layers of two extremal models are perfect conductors and insulators. The responses at the frequency point of intersection are mathematically formulated in the following forms:

(9)\rho (H)={\rho }_a(\omega )[{\displaystyle \frac{\pi }{\mathrm{2}\varphi (\omega )}}-\mathrm{1}],

(10)H=\sqrt{{\displaystyle \frac{{\rho }_a(\omega )}{\omega \mu }}},

Where H is the inversion depth in Eqs. (9) and (10), \rho (H) defines the resistivity, \muis magnetic permeability, \omega is the angular frequency, and {\rho }_a(\omega ) and \varphi (\omega ) are defined as the apparent resistivity and phase at \omega. The frequency-depth in Eq. (10) is similar to the frequency-depth transformation in a half space:

(11)D=\sqrt{{\displaystyle \frac{\rho }{\omega \mu }},}

where D is the exploration depth, \rho is defined as the resistivity of the uniform half space, \mu is the magnetic permeability, and \omega is the angular frequency.

Herein H in Eq. (10) is the depth at the corresponding \omega, and {\rho }_a(\omega ) is considered to be an integrated resistivity up to this depth. To be specific, {\rho }_a(\omega )represents the apparent resistivity at the intersection frequency point of two asymptotic lines. The apparent resistivity at \omega is interpreted based on the impedance calculation. This impedance response is the integrated response covering the whole subsurface below the depth corresponding to\omega. Therefore, in most instances, {\rho }_a in (10) is actually the general resistivity of the subsurface below the corresponding depth of \omega. If the corresponding depth at \omega adequately covered the objective layers, the entire space below this depth could be deemed as a uniform half space with a general resistivity. This general resistivity {\rho }_g, instead of {\rho }_a, should be empirically considered as the compromising resistivity. In addition, (10) is not flexible in terms of thickness extraction of objective layers. To alleviate this problem, we changed the index of the frequency term to interpret the layer thicknesses. The improved equation could be expressed in terms of SI units as follows:

(12)H=\mathrm{356}*\sqrt{{\rho }_g}/f^c,

where H is the estimated depth (m) in (12), {\rho }_g is the general resistivity (W·m), f is the frequency (Hz), and c is the adjustable index of the frequency term. The transformation could be reasonably formulated through trial-and-error tests, suitable for the objective exploration depth and the estimation of the layer thickness. To verify the proposed assessment, we designed three-layer models, and did a comparative analysis between the original and improved frequency-depth transformation in Bostick inversion. To be specific, the models are composed of the 100 W·m cap layer with a thickness of 500 m, an objective layer with the variable resistivity with a 100 m thickness, and a semi-infinite basement layer with the resistivity of 1000 W·m.

For the middle layer of 1 W·m, the actual model and Bostick inversion curves are mapped, as shown in Fig. 5(a). The improved Bostick inversion used a comprising resistivity of 10 W·m covering the objective layer and the subsurface. The indexes were set at 0.1 and 0.3. When the index was 0.1, the thickness identification accuracy was sufficiently high, but the exploration depth encounterd the depth deviation. As for the index of 0.3, the corresponding curve approximated the Bostick inversion curve, but showed an obvious thickness estimation error. Thus, the indexes of the frequency term in the improved frequency-depth transformation could use a number between 0.1 and 0.3, which ascertained that the thickness interpretation accuracy and exploration depth range were both satisfied. The improved Bostick inversion results are depicted in Fig. 5(b) where the middle layer is set to 10 W·m, and the general resistivity to 100 W·m with indexes of 0.3 or 0.5. The improved Bostick inversion results are also shown in Fig. 5(b) illustrating that the high accuracy of thickness identification is achieved when the index is set to 0.3, whereas the depth deviation is revealed. In contrast, another improved result is achieved that is almost identical to the original Bostick inversion, when the index c is set to 0.5. This means that the thickness identification accuracy decreases when the depth reaches the required levels. In this regard, it is suggested that an empirical index between 0.3 and 0.5 be use if the general resistivity is reasonably obtained. It has been proven that if the objective layer is as high as 500 W·m, the high resistive anomaly could not be identified as shown in Fig. 5(c). When the general resistivity is assumed as 600 W·m, and the index is given as 0.5, the improved Bostick inversion could possibly identify an anomalously low valley corresponding to the high resistive layer. This demonstrates that the improved inversion alleviates the identification uncertainty of the original Bostick inversion. However, if the index is set to be 0.7, the high resistive layer is scarcely identified as compared to the original inversion. Thus it should not be assumed that the high resistive anomaly is an interpretation of the inversion.

Practical principle of the SLF technique

The improved Frequency-Depth transformation is reasonably accepted for the SLF technique based on the aforementioned model tests. In a practical geological survey, the subsurface electrical parameters in the study area should be known before carrying out SLF prospecting measurements. The general resistivity and indexes of the frequency terms in the improved Frequency-Depth transformation should be determined based on the model tests or actual geological information.

In order to interpret the objective anomalies, synthetic models were explored by deriving the electrical parameters and SLF response characteristics. Generally speaking, the sedimentary stratigraphy, within a few kilometers of depth, is composed of three parts: sedimentary cap layer, objective layer, and basement layer from top to bottom (Wang et al., 2014a). In Fig. 6(a), the cap layer is shown as a mixture of sandstone and mudstone. Thus, the resistivity could be deemed as 100 W·m. In contrast, the objective layer is observed as more conductive than the surrounding layers according to the aforementioned theoretical analysis, namely, 1 W·m or so. This is commonly observed in geothermal anomalies or coalbed methane reservoirs (Wang et al., 2014b). The basement layer, i.e., the high resistive limestone aquiclude, follows with a resistivity value of 1000 W·m.

As previously demonstrated in the model scales, the cap layer thickness varies from 200 m to 1000 m (i.e., 200 m, 400 m, and 1000 m), and the objective layer is typically 100 m thick over a 1000 W·m uniform basement layer. In the empirical tests, the general resistivity in the Frequency-Depth transformation was set to 400 W·m and the index was set to 0.5. This is reasonably demonstrated in Fig. 6(b). In comparison with the straight trend lines, the upward tails at the corresponding depths of the objective layers are shown by arrows. Namely, if the middle conductive layer is 200 m deep, the upward trend of the tail begins at the depth of 200 m. This demonstrates that the Frequency-Depth transformation is empirically reasonable for the estimated depths of geo-objects. The basement layers are also revealed when the SLF responses become stationary at a depth of 2000 m. The depth ranges of objective layers are inferred to be coarsely located between 200 m and 2000 m.

Therefore, the practical principle of the SLF technique should meet the theoretical requirements. If the objective layers are one or more orders of amplitude more conductive than the surrounding layers, they could be clearly interpreted. When the ratios of “depth to thickness” are shown to be less than eight, the objective depths could then also be interpreted. Moreover, if combined with the synthetic model tests or geological constraints, the Frequency-Depth transformation is empirically obtained, and the depth ranges of geo-objects could be coarsely interpreted.

Applicability analysis of the SLF technique

It is essential to acquire a high-quality magnetic amplitude spectrum over a broad frequency range given the benefit derived from the high-quality data, such as improved depth resolution and target imaging. The non-imaging SLF remote sensing system is composed of three parts: 1) a high-resolution remote sensor for sensing the magnetic field, 2) a processing unit for controlling acquisition, digitization, filtering, and recording of raw data, and 3) the power supply. This instrumentation has been developed and continuously improved for sufficient sensitivity and extensive ranges for measuring natural source broadband SLF signals. The sensor shown in Fig. 7 (a) is a small,portable induction coil with low power consumption. This type of sensor is required due to its sensitivity to low-frequency signals, which can be increased by inserting a pre-amplification circuit. The signals can then be further processed by multi-scale amplifying-filtering, A/D transformation, and Frequency-Depth transformation modules controlled bythe processing unit shown in Fig. 7 (b). The “Depth-Magnetic Relative Amplitude” curves are ultimately displayed in real-time on the screen of the processing unit. The power consumption of the remote sensing system is less than 15 w, and as such, a portable power supply (Fig. 7 (c)) is needed to meet the requirements of conducting field surveys. The stability of the SLF instrumentation has been verified by reliable experimental and observational tests and has been shown to be effective and economical for conducting regional geological surveys.

A detailed working flowchart of the SLF technique is summarized and illustrated in Fig. 8. The practice is composed of a ground survey, data acquisition, data processing, and geological interpretation. The appropriate sites in a specific study area should be deployed according to the ground survey. To reduce unnecessary interference, the ground survey is of utmost importance for site deployment. SLF sites should be relatively flat to enable accurate positioning of the SLF system and to minimize the noise level due to ground vibration. Therefore, artificial noise sources, such as power lines and vehicles, should be at a significant distance from the SLF sites. In addition, the data acquisition process should meet the requirements for minimizing motion-induced noises. The acquisition parameters, such as detector direction, step interval, scale amplification, and depth ranges, should be defined to ensure that the results optimally follow project aims. The SLF system is then utilized to repetitively acquire natural source signals in a multi-angle manner. To ensure that the response errorbars are within 10%, the SLF responses of each angle are considered throughout the correlation analysis and repeatability evaluation (Wang et al., 2014a).

An SLF processing module has been developed for pre-processing, quality evaluation, and de-noising. The module contains functions of format transformation, stacking, noise filtering, and normalization. All the raw data would be pre-processed and evaluated to obtain the high SNR (signal-noise-ratio) and repetitive SLF magnetic amplitude spectrum responses with adequately high correlations. Furthermore, the SLF curves should be normalized and processed by de-noising algorithms to preliminarily filter both random and artificial noise (cultural noise). Data processing methods suitable for this type of non-stationary SLF signals are also suggested, such as Independent Component Analysis (ICA), Wavelet and Lifting Wavelet analysis, as well as Hilbert-Huang Transform (HHT) (Wang et al., 2014a). SLF curves are transformed to the depth domain using the specific Frequency-Depth transformation. The ultimately processed SLF curve is then ready for geological interpretation of objective layers. In a word, the applicability of the SLF technique is not only theoretically satified, but also determined by convincing data acquisition, processing, and interpretation.

Case studies

Thermal reservoir prospecting

Geothermal energy, which can effectively reduce the world’s dependence on fossilized energy in favor of environmental protection (Meju, 2002), is one of the most attractive energy resource options across the globe. Geothermal reservoirs normally dissolve massive mineralized ions at high temperatures, and are more conductive than the adjacent layers. Generally, the resistivity is one order of magnitude lower than the background, so that the strong resistive discrepancy provides the basis for reservoir identification using the SLF technique (Heise et al., 2008; Spichak and Manzella, 2009).

Geological settings

The geothermal area at Peking University, which is currently under investigation, lies at the northeastern Mentougou fold trap of the Yanshan platformal fold belt, close to the core in the northern limb of the Kunming Lake anticline. Figure 9(a) depicts the entire geothermal area where the Huangzhuang-Gaoliying fault and Babaoshan fault are located in the southeast. This is a regional extension fault, and is the primary structure suited for heat and water migration. Babaoshan fault is a reverse fault and increases the burial depths of reservoirs. Both faults account for the favorable areas for thermal resource accumulations where geothermal wells are widely distributed throughout the central and the northern sections (shown by black points).

According to the drilling data depicted in Fig. 9(b), the main geothermal reservoir is in the Wumishan Form, Jixian System (J_xw). The overlying layers are include the Quaternary system (Q), Ordovician system (O), Cambrian system (∈), Qingbaikou System (Q_n), and Jixian System (J_x) in a sequence. The depth range of 1.6 km in Q, O, and ∈ systems is the thick sedimentary cap layer, where the average temperature is about 27°C. It has a large thermal contrast with the underlying Q_n system with a temperature of 41°C. The Q_n system seals the geothermal regions below, and has almost the same lithology as those of geothermal reservoirs. The Tieling Form (J_xt) and Hongshuizhuang Form (J_xh) in the Jixian System (J_x) are the principal cap layers for reservoirs, and meet the requirements of geothermal reserves. The entire Jx system is composed of dolomite, shale, slate, and quartzite with a maximum average temperature of 60°C. The geometrics of the study area is geologically plausible as a 1-D hypothetical model, including the thick cap layer, the objective geothermal reservoir, and the dolomite basement. This case can be used for investigating the application of the SLF technique.

Reservoir prospecting of Well jr168

Well jr168 is located in Chengfu Garden, Peking University, shown in Fig. 10. The outlet temperature reaches a maximum of 60°C, and the water yield is 2000 m³ per day. The total thickness of the cap layer over the geothermal reservoir (J_xw Form) is 2630 m covering many water-rich segments (Qin et al., 2010). The depth range of the thermal reservoir is 2630.00 ‒3218.68 m, and the reservoir temperature is 56.6°C‒65.60°C. The drilling information is shown in Table 1. This infers that the geothermal reservoir is a low resistive layer, and the ratio of “depth to thickness” is nearly 4. According to the model tests discussed in Section 3, it meets the requirements of reservoir identification.

The SLF prospecting survey was carried out in November 2013, and the acquisition site was selected 30‒50 m far away from each artificial source in order to reduce the cultural noise. The SLF prospecting system was switched every other 90 degrees, and acquired SLF data five times for multiple angles. The collected data were evaluated for quality assurance and were normalized. The two empirical parameters in the Frequency-Depth transformation were set to 1000 W·m and 0.5. A raw SLF curve (illustrated in Fig. 10(a))was later processed by the empirical mode decomposition (EMD), curve reconstruction, and wavelet de-noising (Wang et al., 2014a). As shown in Fig. 10(b), the ultimate SLF curve was used for the following reservoir interpretation with 50 Hz powerline noise and random noise suppressed.

Reservoir interpretation result

The reservoir interpretation result,shown in Fig. 11, demonstrates that the SLF curves could be fitted with a cubic polynomial (Root-Mean-Square (RMS)>90%) in a satisfactory agreement (Wang et al., 2014a) . This fitting of a curve reveals the main features of the original SLF curve. We found that the curve changed slightly and there was no conductive layer above the depth of 1700 m (shown by the first arrow), indicating the turning point of the fitting curve as compared to the trend line. The interpreted depth range between 1000 m and 1700 m could be related to that of the Cambrian System (ϵ), as shown in Table 1. The logging data showed it as a resistive layer (300–2000 W·m) that is impregnated with the main component of carbonate rocks. This is also explained by a slightly changing underground temperature, with an average warming rate of only 0.63°C/100 m. In contrast, there is a small uprising amplitude between 1700 m and 2200 m, shown by the first and second black arrows. This is correspondingly interpreted as the conductive Qingbaikou System (Q_n). The drilling data in Table 1 reveal that the average warming rate can reach up to 2.45°C/100 m. Meanwhile, the logging data show an average resistivity of 50–100 W·m for the Q_n system, which can be explained by the active electrons and aqueous ions in the hot environment where the layer is more conductive. Under the depth of 2200 m (shown by the third arrow), there is a larger uprising amplitude in the fitting curve compared to the trend line (the dotted line). The lower resistive distribution between 2200 m and 3200 m is inferred as the Jixian System (J_x) according to the drilling data. Given that the tempreture rapidly rises at the rate of 2.01°C/100 m from 2279 m to 3218 m, the average resistivity of this system revealed by logging data is lower than 30 W·m. Interpreted geothermal distribution is in perfect consistency with the anomalous conductive area between 2200 m and 3200 m. Although the actual depth range of the commercial geothermal reservoirs is 2630.00–3218.68 m in Jixian system, the aquifers, as adjacent layers, are more widely spread. However, the general distribution of geothermal reservoirs was still determined by imaging the conductive area. There might be a high resistive basement below the depth point of 3200 m, which is the interaction point of the trend line and the fitting line. The subsurface space is inferred as the high resistive basement layer owing to the gradual curve of the SLF interpretation which approaches to a stable value from the depth of 3218 m. This is one of the observed consistent outcomes of known model tests. This study demonstrates that the SLF technique is a promising tool for directly interpreting conductive thermal reservoirs.

Coalbed methane reservoir interpretation

Coalbed methane (CBM), a complimentary resource and trigger for a revolution in traditional energy structures (Moore, 2012), is one of the most significant, unconventional energy resources worldwide. The SLF technique has preliminarily been applied in the identification and monitoring of CBM reservoirs (Wang et al., 2014a). As an example, we chose a study area in Qinshui Basin, China, which was one of the most significant commercial CBM fields with a total reserve of 75.4×10⁹ m³ (Su et al., 2005).

Study area

The study area is located in the southeastern section of Qinshui Basin, dipping 5°NW. The geological setting within the area is relatively simple and flat, plausible as a 1-D geo-electrical structure. Few faults exist, but folds with a striking direction of NNE-SSW or N-S are common. According to the previous investigations, the sedimentary layers, in sequence, include the Quaternary system (Q), Shangshihezi Form, Permian system (P₂s), Xiashihezi Form, Permian system (P₁x), Shanxi Form, Permian System (P₁s), Taiyuan Form, Carboniferous system (C₃t) and basement layers as depicted in Fig. 12(a) (Liu et al., 2014). The main CBM units of Shanxi and Taiyuan Forms are over 150 m thick, and the economic coal seams are No. 3 and No. 15. The No. 3 coal seam is located in Shanxi Form, P₁s, and No. 15 coal seam is located in Taiyuan Form, C₃t. The total coal-bed thickness is over 15 m. The bottom depth of the No. 15 coal seam is less than 650 m, and the No. 3 coal seam is 20 m to 50 m shallower than the top depth of No. 15 (Su et al., 2005).

Figure 12(b) shows that the gas content of the No. 3 coal seam is generally higher than 10 m³/t. All of the reservoirs are hydraulically blocked by low resistive aquifers, with the top depth starting at 460‒500 m. According to the ground survey, the CBM wells are pumping massive groundwater from the reservoirs while producing CBM. This means the CBM reservoirs are filled with groundwater, and the aquifers, including the reservoirs, form an entire conductive area. Therefore, the total thickness of aquifers is inferred as 100‒200 m. According to the drilling data, the resistivity of the aquifers, including CBM reservoirs, is lower than 50 W·m following the water-gas production. Thus, the ratio of “depth to thickness” is 2.5–5, and reaches the depth resolution of conductive anomalies. In addition, non-producing adjacent layers are normally 200–1000 W·m. On the basis of the aforementioned feasibility analysis, configuration discrepencies between reserovirs and the background contribute to an indirect resolution for low resistive aquifers using the SLF technique.

We carried out the SLF prospecting measurements around Well 01 and Well 02 in May 2010. In order to reduce the noise level, the prospecting systems were 50 m away from the production wells and other artificial sources. The direction of the sensor was rotated every other 90 degrees at each site, and the data acquisition was performed five times for each angle. To reduce the response errorbars within 10%, the SLF responses of each angle were continually evaluated throughout. All the SLF data were processed and transformed according to the working flowchart (Fig. 8). Noise was reasonably suppressed using empirical mode decomposition and wavelet de-noising methods (Wang et al., 2014a). As observed through field survey and model tests, the general resistivity of the Frequency-Depth transformationwas 400 W·m and the index was 0.5 in this study area. The ultimately processed SLF curves were employed for reservoir interpretation.

Results and discussions

The SLF interpretation curve of Well 01 is depicted in Fig. 13 (a). It is inferred that the turning point was at a depth of 480 m, where the curve had an obvious uprising amplitude in comparison to the trend line (dotted line). This demonstrates that the top depth of the conductive layer was detected. Considering the exploration depth range, a potential conductive anomaly was interpreted above 600 m. The thickness of aquifers was inferred over 100 m. According to the drilling data of Well 01, No. 15 coal seam is 545‒551 m. Thus, the interpreted conductive anomalies were associated with No. 15 reservoir and the adjacent aquifers. The depth range of No. 15 reservoir was not identified from the uprising amplitude of the SLF interpretation curve. However, it is worth noting that there was an obvious high amplitude wave packet in the SLF curve (pointed out by the dotted box). With the artificial and random noise suppressed, the anomalous packet could be low-frequency EM radiation signals from producing CBM reservoirs (Wang et al., 2014b). According to the previous research, EMR signals were generated in gas-bearing coal deformation or fracturing. Different types of coal samples would generate EM radiation signals when CBM went through them (He et al., 2012). Experiments demonstrate that if the gas pressure is increased, consequently, the radiation intensity also increases (Liu and He, 2001; Frid and Vozoff, 2005). EMR intensity could be positively enhanced with increased fracture strength, and quantitative coupling correlations could also be implemented between coal or rock uniaxial load and EMR amplitudes (He et al., 2012). EMR signals commonly distribute at a wide frequency range (Wang et al., 2014), with diverse signal frequencies being motivated at low frequency (Wang et al., 2011). The signals in the frequency range of 10 kHz or higher, are relatively too weak to be observed on the Earth’s surface. Due to the “skin depth” theory, high-frequency EMR signals travel a few hundred meters below the surface, and possibly attenuate to zero when they are acquired. Thus, EMR signals at frequencies lower than kHz would be primarily reserved for revealing information on subsurface anomalies. Wang et al. (2014) investigated the ultra-low-frequency (300‒3000 Hz, i.e., ULF) EMR signals in coal mines to monitor reservoir stress and failure (Wang et al., 2014). In spite of the coexistence of CBM and water, low-frequency EMR signals are also generated when CBM flows through the coal samples (Frid et al., 2003). Advanced data processing methods to enhance the low frequency (3‒3000 Hz) signals were developed, and the EMR over or greater than 1 kHz were further suppressed (Wang et al., 2013). The EM radiation signals diffuse in pulses which stack to form a local wave packet in the frequency spectrum. If the Frequency-Depth transformation is used, the wave packets could correspond to the anomalies at a specific depth (Wang et al., 2014b). In addition, the intensity of radiation signals observed on the surface normally exceeds the magnitude of mV, which could be observed by the high-resolution prospecting system (Greiling and Obermeyer, 2010). If low-frequency radiation signals from the producing reservoirs are exceptionally intensive, the frequency-corresponding depth distributions of reservoirs could be interpreted. In Fig. 13 (a), the anomalous high-amplitude packet is interpreted at the depth of 540‒550 m. It is located in correspondence to the position of No. 15 coal seam with a depth range of 545‒551 m. Additionally, Well 01 was indeed in the rapid production status, which was inferred to account for the origin of radiation signals. Therefore, the anomalous radiation packet extracted from the SLF curve is consistent with the actual production circumstances.

Well 02 is located 100‒200 m the east of Well 01. Geological information reveals that fault F1 is located between the two wells. This fault explains the dropdown of aquifers and reservoirs, and thus, the low resistive area which starts at a depth of 510 m. The No. 15 reservoir, which has been observed at a depth of 560–570 m, is 20 m deeper than Well 01. This indicates that the fault displacement is about 20–30 m. In contrast to the SLF response (Fig. 13 (a)), the uprising amplitude in the deep SLF response of Well 02 is less than that of Well 01 (Fig. 13(b)). This shows a higher resistivity for No. 15 reservoir than for Well 01. According to the production status of the two wells, the gas production rate of Well 01 was obviously slower than Well 02. This indicates that No. 15 reservoir of Well 01 accumulated more groundwater, and produced more water than gas.As a result, gas production gradually decreased until entering the stopping period. In contrast, the water content in the No. 15 reservoir of Well 02 was relatively lower resulting in a rapid increase in gas production. These results are consistent with the more intensive radiation signals (i.e., high amplitude packets in dotted box) as shown in Fig. 13(b). This study demonstrates that the SLF technique provides an attractive tool in reservoir identification and dynamic monitoring.

Conclusions

In this paper, the Super-Low Frequency EM prospecting technique is proposed as a promising geological survey tool. The theoretical basis, numerical algorithms, feasibility analysis, practical principle, and geological application have all been illustrated.

1) The theoretical MT and SLF responses in 1-D, 2-D, and 3-D models in the frequency domain were consistent, which verifies that our numerical algorithms are reliable and convincing.

2) The modeling results of real 3-D models show that averagely normalized SLF magnetic amplitudes should be employed in geological interpretation, which also suggests that multi-angle data acquisition should be performed many times for each site in practical application.

3) Three-layer model results demonstrate that the SLF technique is more sensitive to the conductive objective anomalies, especially when the resistivity of objective layers is one order or more of amplitude lower than adjacent layers. The resistive layers are reluctant to be interpreted.

4) An improved Frequency-Depth transformation was proposed in order to semi-quantitatively accommodate the objective depth and objective thickness by adjusting the general resistivity and frequency term indexes. Consequently, a complete workflow was proposed for practical geological survey.

5) The case studies in geothermal exploration and coalbed methane reservoir interpretation show that the SLF technique is suitable for imaging the distributions of conductive anomalies, and meanwhile, it is in favor of producing reservoir identification if the EM radiation packets are extracted. Furthermore, the production status of reservoirs is also consistent with the SLF interpretation results.

This study lays a solid foundation for the promising feasibility of the SLF technique in subsurface geological exploration.

References

[1] Aprea C, Booker J, Smith J (1997). The forward problem of electromagnetic induction: accurate finite-difference approximations for two-dimensional discrete boundaries with arbitrary geometry. Geophys J Int, 129(1): 29–40

[2] Atzberger C (2013). Advances in remote sensing of agriculture: context description, existing operational monitoring systems and major information needs. Remote Sens, 5(2): 949–981

[3] Avdeev D B (2005). Three-dimensional electromagnetic modelling and inversion from theory to application. Surv Geophys, 26(6): 767–799

[4] Berdichevsky M N, Dmitriev V I, Zhdanov M S (2010). Possibilities and problems of modern magnetotellurics. Izv Phys Solid Earth, 46(8): 648–654 doi:10.1134/S1069351310080021

[5] Börner R U (2010). Numerical modelling in geo-electromagnetics: advances and challenges. Surv Geophys, 31(2): 225–245 doi:10.1007/s10712-009-9087-x

[6] Congalton R, Gu J, Yadav K, Thenkabail P, Ozdogan M (2014). Global land cover mapping: a review and uncertainty analysis. Remote Sens, 6(12): 12070–12093 doi:10.3390/rs61212070

[7] Cui R B, Qin Q M, Li B S, Wang Q P (2009). Design and development of passive super low frequency electromagnetic data processing software. In: 2009 IEEE International Geoscience And Remote Sensing Symposium, Vols 1-5 IEEE International Symposium on Geoscience and Remote Sensing IGARSS. (345 E 47TH ST, New York, NY 10017 USA: IEEE), 627–630

[8] Frid V, Rabinovitch A, Bahat D (2003). Fracture induced electromagnetic radiation. J Phys D Appl Phys, 36(13): 1620–1628

[9] Frid V, Vozoff K (2005). Electromagnetic radiation induced by mining rock failure. Int J Coal Geol, 64(1‒2): 57–65

[10] Gomez-Treviño E (1996). Approximate depth averages of electrical conductivity from surface magnetotelluric data. Geophys J Int, 127(3): 762–772 doi:10.1111/j.1365-246X.1996.tb04055.x

[11] Greiling R O, Obermeyer H (2010). Natural electromagnetic radiation (EMR) and its application in structural geology and neotectonics. J Geol Soc India, 75(1): 278–288

[12] He X Q, Nie B S, Chen W X, Wang E Y, Dou L M, Wang Y H, Liu M J, Hani M (2012). Research progress on electromagnetic radiation in gas-containing coal and rock fracture and its applications. Saf Sci, 50(4): 728–735

[13] Heise W, Caldwell T G, Bibby H M, Bannister S C (2008). Three-dimensional modelling of magnetotelluric data from the Rotokawa geothermal field, Taupo Volcanic Zone, New Zealand. Geophys J Int, 173(2): 740–750

[14] Huang H P, Won I J (2003). Detecting metal objects in magnetic environments using a broadband electromagnetic method. Geophysics, 68(6): 1877–1887 doi:10.1190/1.1635040

[15] Liu H, Sang S, Wang G, Li M, Xu H, Liu S, Li J, Ren B, Zhao Z, Xie Y (2014). Block scale investigation on gas content of coalbed methane reservoirs in southern Qinshui basin with statistical model and visual map. J Petrol Sci Eng, 114: 1–14

[16] Liu M J, He X Q (2001). Electromagnetic response of outburst-prone coal. Int J Coal Geol, 45(2-3): 155–162

[17] Mackie R L, Madden T R, Wannamaker P E (1993). Three-dimensional magnetotelluric modeling using difference equations—Theory and comparisons to integral equation solutions. Geophysics, 58(2): 215–226

[18] Meju M A (2002). Geoelectromagnetic exploration for natural resources: Models, case studies and challenges. Surv Geophys, 23(2‒3): 133–206 doi:10.1023/A:1015052419222

[19] Melesse A M, Weng Q, Thenkabail P S, Senay G B (2007). Remote sensing sensors and applications in environmental resources mapping and modelling. Sensors (Basel Switzerland), 7(12): 3209–3241

[20] Moore T A (2012). Coalbed methane: a review. Int J Coal Geol, 101: 36–81

[21] Ouchi K (2013). Recent trend and advance of synthetic aperture radar with selected topics. Remote Sens, 5(2): 716–807

[22] Qin  Q M, Li  B S, Cui  R B, Jiang  H B, Wang  Q P, Zhang  Z X (2010). Analysis of factors affecting natural source SLF electromagnetic exploration at geothermal wells. Chinese Journal of Geophysics, 53(3): 685–694doi: 10.3969/j.issn.0001-5733.2010.03.023(in Chinese)

[23] Rodríguez J, Esparza F J, Gómez-Treviño E (2010). 2-D Niblett-Bostick magnetotelluric inversion. Geologica Acta, 8(1): 15–30 doi: 10.1344/105.000001513

[24] Schmucker U (1987). Substitute conductors for electromagnetic response estimates. Pure Appl Geophys, 125(2‒3): 341–367

[25] Shan Z, Qi Q, Nan W, Li C, Yan B, Cheng Y (2014). Analysis and design of passive super low frequency detection system. In: 2014 IEEE Geoscience and Remote Sensing Symposium. (IGARSS). Proceedings (Piscataway, NJ, USA: IEEE), 5056–5059

[26] Siripunvaraporn W, Egbert G, Lenbury Y (2002). Numerical accuracy of magnetotelluric modeling: a comparison of finite difference approximations. Earth Planets Space, 54(6): 721–725 doi:10.1186/BF03351724

[27] Slonecker T, Fisher G B, Aiello D P, Haack B (2010). Visible and infrared remote imaging of hazardous waste: a review. Remote Sens, 2(11): 2474–2508

[28] Spichak V, Manzella A (2009). Electromagnetic sounding of geothermal zones. J Appl Geophys, 68(4): 459–478

[29] Su X B, Lin X, Liu S, Zhao M, Song Y (2005). Geology of coalbed methane reservoirs in the Southeast Qinshui Basin of China. Int J Coal Geol, 62(4): 197–210

[30] van Genderen J L (2011). Advances in environmental remote sensing: sensors, algorithms, and applications. Int J Digit Earth, 4(5): 446–447

[31] Wang E Y, He X Q, Liu X, Li Z, Wang C, Xiao D (2011). A non-contact mine pressure evaluation method by electromagnetic radiation. J Appl Geophys, 75(2): 338–344

[32] Wang E Y, Jia H L, Song D Z, Li N, Qian W H (2014). Use of ultra-low-frequency electromagnetic emission to monitor stress and failure in coal mines. Int J Rock Mech Min Sci, 70: 16–25

[33] Wang N, Qin Q M, Chen L, Bai Y B, Zhao S S, Zhang C Y (2014a). Dynamic monitoring of coalbed methane reservoirs using Super-Low Frequency electromagnetic prospecting. Int J Coal Geol, 127: 24–41

[34] Wang N, Qin Q M, Chen L, Bai Y B, Zhao S S, Zhang C Y (2014b). Passive Super-Low Frequency Remote Sensing Technique For Monitoring Coal-Bed Methane Reservoirs. In: 2014 IEEE International Geoscience And Remote Sensing Symposium (IGARSS) IEEE International Symposium on Geoscience and Remote Sensing IGARSS. (345 E 47TH ST, NEW YORK, NY 10017 USA: IEEE), 867–870

[35] Wang N, Qin Q M, Xie C, Chen L, Bai Y B (2013). Coal-Bed Methane Reservoir Identification Using The Natural Source Super-Low Frequency Remote Sensing. In: 2013 IEEE International Geoscience And Remote Sensing Symposium (IGARSS) IEEE International Symposium on Geoscience and Remote Sensing IGARSS. (345 E 47TH ST, NEW YORK, NY 10017 USA: IEEE), 4022–4025

[36] Yavuz M E, Teixeira F L (2009). Ultrawideband microwave sensing and imaging using time-reversal techniques: a review. Remote Sens, 1(3): 466–495

[37] Zhdanov M S, Varentsov I M, Weaver J T, Golubev N G, Krylov V A (1997). Methods for modelling electromagnetic fields Results from COMMEMI—the international project on the comparison of modelling methods for electromagnetic induction. J Appl Geophys, 37(3‒4): 133–271 doi:10.1016/S0926-9851(97)00013-X

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg