TY - JOUR
AU - Karimi,, Kurosh
AB - Abstract Unquestionably, current density pseudo-section derived from in-phase very low frequency (VLF) data, developed by Karous and Hjelt, has been one of the most popular preliminary schemes for the interpretation of VLF data due to its straightforwardness and applicability. In practice, however, the Karous and Hjelt filtering method does not behave accurately and has some weaknesses: (i) the depth estimation is erratic, even in shallow cases, (ii) to a far lesser extent, the horizontal estimation is inaccurate, in particular when the sources interfere with each other, and (iii) under circumstances where multiple bodies are in close proximity, spurious artifacts appear. These factors could be misleading and result in a wrong interpretation. Changing the Biot–Savart law into a linear equation and employing inverse theory in this paper, the precision of spatial estimations, in the case of superficial objects, becomes far more precise. More importantly, independent of the number of neighboring targets and their closeness, no spurious artifacts are observed. This methodology is tested on synthetic isolated and interfering various models with different sizes, depths, shapes and resistivities. It is also applied to real data over the ‘Ghoori Ghale’ cave in the vicinity of Kermanshah, Iran, and its satisfying and noteworthy solutions are compared with the Karous and Hjelt method. VLF, Karous–Hjelt filtering, least squares inversion, current density 1. Introduction The very low frequency (VLF) method is an electromagnetic induction-based method that has been used in a variety of purposes like mapping of mineralization in fault zones (Philips & Richards 1975), exploration of mineralized bodies (Parker 1980), nuclear waste detection (Sinha & Hayles 1988; Ramesh Babu et al.2004), ground water detection or contamination (Monterio Santos et al.2006; Sundararajan et al.2007), soil engineering and detection of basement fractures (Ramesh Babu et al.2007), archaeological investigations (Khalil et al.2010; Timur 2012) etc. Two conventional types of instruments, i.e. Geonics EM-16 (Paterson & Ronkav 1971) and SDVR-3 (Sedelnikov et al.1971), measure the dip of ellipse major axis and VLF field components, respectively. In the latter case, the ratio of the vertical anomalous field, Hz, to the horizontal primary field, H0, has been elaborately studied and it has been shown that in the presence of conductors or conducting zones, the in-phase component of this ratio indicates crossovers (McNeill & Labson 1991). There is much early evidence of the utility and practicality of the VLF technique in the literature (e.g. Gharibi & Pedersen 1999; Khalil & Monterio Santos 2010; Sundararajan et al.2011; Al Bulushi et al.2016). Nonetheless, this survey focuses on the current density pseudo-section. Karous & Hjelt (1977, 1983) (K&H) proposed a filtering method and plotted a current density pseudo-section map to specify a rough location of shallow anomalous sources. However, their technique had two main problems: (i) when dealing with relatively small targets, even in shallow positions, high current spots barely coincided with the true location of conductive bodies and (ii) in the presence of multiple sources, spurious artifacts emerged (Ogilvy & Lee 1991). The K&H filter was constructed based on the magnetic fields caused by currents flowing beneath the surface of the ground, and was more comprehensive than the filter developed by Fraser (1969). Sundararajan et al. (2006) compared and contrasted the Fraser and K&H filters in detail. Four years later, Gnaneshwar et al. (2011) employed an in-phase VLF anomaly component and studied geological boundaries and shear zones of the Schirmacheroasen area in east Antarctica by means of an analytic signal approach as well as Fraser and K&H filterings. In the same year, the difference between the responses of K&H filtering and Inv2DVLF software was evaluated (Khalil & Monterio Santos 2010). In that survey, Khalil & Monterio Santos (2010) showed that both techniques yielded very similar results, except that the K&H technique had a depth limitation. In this study, the anomalous bodies are assumed to be 2D and elongated perpendicular to the profile directions. The discrepancy between solutions of K&H filtering and the inverse method is perused. It turns out that the inverse method is more trustworthy than the filtering one, specifically in the case of superficial anomalies. 2. Mathematical tools Here, we modify the idea presented by Karous & Hjelt (1983) for the interpretation of VLF data. Our method uses the least squares inversion scheme to recover a current density section beneath a profile over which the data are acquired. Consider a current density J, which is induced by a uniform magnetic field, H0: the current flows in an element of a subsurface continuous conductive body spanning horizontally from a to b and vertically from Z1 to Z2 (figure 1). The ratio of the vertical component of the secondary magnetic field due to induced currents, Hz, to the primary magnetic field, H0, can be simply calculated by the Biot–Savart law as (modified from Reitz & Milford 1966): $$\begin{equation}\frac{{{H_z}}}{{{H_0}}} = \frac{1}{{2\pi {H_0}}}\int_{a}^{b}\!{{\int_{{{Z_1}}}^{{{Z_2}}}{{\frac{{(x - t)}}{{[{{(x - t)}^2} + {Z^2}]}}\,\,J(t)\,\,dz}}\,\,dt}}.\end{equation}$$ (1) Figure 1. Open in new tabDownload slide Schematic plan of a ground cross-section along with unit cells bearing electrical currents, depth layers and data points. I’q is the current in the qth depth layer and I’(q-1)n+s is the electrical current related to the ‘(q − 1)n + s’th unit cell, where s is between 1 and n. Figure 1. Open in new tabDownload slide Schematic plan of a ground cross-section along with unit cells bearing electrical currents, depth layers and data points. I’q is the current in the qth depth layer and I’(q-1)n+s is the electrical current related to the ‘(q − 1)n + s’th unit cell, where s is between 1 and n. In equation (1), x is the data point coordinate, t is the horizontal location of a point inside the conductive element and J(t) is the true current density. According to equation (1), the current density J is not a function of x or z, i.e. the equation is explicit versus J. So, one can easily separate J from the rest of the integrand. Also, the equation (1) is equal to the convolution of J and x/(x2 + z2) and then, because of the explicit role of J, the integral convolution can be transformed into a discreet convolution. Let us denote Hz/H0 with just an H. Considering |$g(x,t) = \frac{1}{{2\pi }}\int_{{{Z_1}}}^{{{Z_2}}}{{\frac{{(x - t)}}{{[{{(x - t)}^2} + {Z^2}]}}dz}}$| and |$I(t) = \frac{{J(t)}}{{{H_0}}}$|⁠, One can write: $$\begin{equation}H(x) = \int_{a}^{b}{{g(x,t)\,I(t)\,dt}}.\end{equation}$$ (2) In a discrete form of measurement, the obtained H in the ith measurement point could be indicated as: $$\begin{equation}{H_i} = \int_{a}^{b}{{{g_i}(t)\,I(t)\,dt}},\quad i = 1,2,\ldots,n.\end{equation}$$ (3) Likewise, if the subsurface is divided into separate blocks with discrete center coordinates (tj) and dimensions of ΔZ and ΔT (we suppose that ΔT = ΔZ), which represent the depth extent and width of each block, respectively, equation (2) turns into: $$\begin{eqnarray}{H_i} &&= \sum\limits_{j = 1}^m {[{g_i}({t_j})\Delta T]\,I({t_j})}\\ &&= \sum\limits_{j = 1}^m {{G_i}({t_j})\,I({t_j})},\left\{ {\begin{array}{@{}*{1}{l}@{}} {i = 1,2,...,n}\nonumber\\ {j = 1,2,...,m} \end{array}} \right.\end{eqnarray}$$ (4) where $$\begin{eqnarray*}{G_i}({t_j}) &=& \left( {\frac{{\Delta Z}}{{2\pi }}\frac{{({x_i} - {t_j})}}{{[{{({x_i} - {t_j})}^2} + {{(Z)}^2}]}}}\right)\\ &&\Delta T = \left( {\frac{{\Delta Z}}{{2\pi }}\frac{{[(i - 1) - (j - 1)]}}{{{{[(i - 1) - (j - 1)]}^2} + {{(q)}^2}}}}\right),\end{eqnarray*}$$|${x_i} = a + (i - 1)\Delta T, {t_j} = a + (j - 1)\Delta T$| and |$\Delta T = \frac{{b - a}}{n}$|⁠. Notice that the depth of each row of the current blocks (depth layer) is Z = q ΔT, in which q = 1, 2, … L (figure 1). Up to now, the formulations have been similar, to some extent, to those from Karous and Hjelt (1983). In a system with n measured data, equation (4) can be written as $$\begin{equation}{\boldsymbol{H}} = {\boldsymbol{GI}}.\end{equation}$$ (5) Here, G is an n by m matrix, where n is the number of columns (refer to figure 1) and m = (n − 1)L and L is the number of depth rows (depth layers). As is clear, in the discrete condition, equation (1) converts into a linear relation. Equation (5) resembles the well-known fundamental equation of the inverse theory as (Aster et al.2003): $$\begin{equation}{\boldsymbol{d}} = {\boldsymbol{Gm}}.\end{equation}$$ (6) In this case, we deal with an underdetermined system of linear equations. The damped least squares solution to extract I is (modified from Tikhonov & Arsenin 1977): $$\begin{equation}{{\bf I = [}}{{{\bf G}}^{{\bf T}}}{{\bf G}} + \beta {I_d}{]^{ - 1}}{{{\bf G}}^{{\bf T}}}{{\bf H}},\end{equation}$$ (7) where Id is an identity matrix and β is a positive damping factor. In fact, β is a trade-off parameter between minimization of data and model objective functions (Aster et al.2003). In equation (7), |${{\bf G = [G}}_1^{\prime}\,\,{{\bf G}}_2^{\prime}\,\,{{\bf G}}_3^{\prime}\,\,...\,\,{{\bf G}}_q^{\prime}\,\,...\,\,{{\bf G}}_L^{\prime}{{\bf ]}}$|⁠, |${{\bf H = }}{[ {{{{\bf H}}_{{\bf 1}}}\,\,{{{\bf H}}_2}\,\,{{{\bf H}}_3}\,\,...\,\,{{{\bf H}}_n}} ]^T}$| and |${{\bf I = }}{[ {{{\bf I}}_1^{\prime}\,\,{{\bf I}}_2^{\prime}\,\,...\,\,{{\bf I}}_q^{\prime}\,\,...\,\,{{\bf I}}_L^{\prime}} ]^T}$| where T denotes transpose and $$\begin{eqnarray*} {{\bf G}}_{{\bf q}}^{{\bf ^{\prime}}} = \left[ {\begin{array}{@{}*{5}{c}@{}} {G_{1,(q - 1)n + 1}^{\prime}}&\,\,{G_{1,(q - 1)n + 2}^{\prime}}&\,\,{G_{1,(q - 1)n + 3}^{\prime}}&\,\,{...}&{G_{1,nq}^{\prime}}\\ {G_{2,(q - 1)n + 1}^{\prime}}&\,\,{G_{2,(q - 1)n + 2}^{\prime}}&\,\,{G_{2,(q - 1)n + 3}^{\prime}}&\,\,{...}&{G_{2,nq}^{\prime}}\\ {G_{3,(q - 1)n + 1}^{\prime}}&{G_{2,(q - 1)n + 2}^{\prime}}&{G_{3,(q - 1)n + 3}^{\prime}}&\,\,{...}&{G_{3,nq}^{\prime}}\\ \begin{array}{@{}l@{}} :\\ . \end{array}&\,\,\begin{array}{@{}l@{}} :\\ . \end{array}&\,\,\begin{array}{@{}l@{}} :\\ . \end{array}&\,\,\begin{array}{@{}l@{}} :\\ . \end{array}&\,\,\begin{array}{@{}l@{}} :\\ . \end{array}\\ {G_{n,(q - 1)n + 1}^{\prime}}&\,\,{G_{n,(q - 1)n + 2}^{\prime}}&\,\,{G_{n,(q - 1)n + 3}^{\prime}}&\,\,{...}&\,\,{G_{n,nq}^{\prime}} \end{array}} \right] \quad {\rm and}\quad{{\bf I}}_{{\bf q}}^{{\bf ^{\prime}}}{{\bf }}=\left[ {\begin{array}{@{}*{1}{l}@{}} {{{{{\bf I^{\prime\prime}}}}_{(q - 1)n + 1}}}\\ {{{{{\bf I^{\prime\prime}}}}_{(q - 1)n + 2}}}\\ {{{{{\bf I^{\prime\prime}}}}_{(q - 1)n + 3}}}\\ :\\ :\\ {{{{{\bf I^{\prime\prime}}}}_{qn}}} \end{array}} \right].\end{eqnarray*}$$ 3. Synthetic models In order to test the inverse method and evaluate its accuracy compared with the K&H filter, several isolated and multiple models with different sizes, depths, shapes and resistivities, with and without overburden in host rocks with varying resistivities, are studied. The models are isolated dipping dike, isolated synclinal structure, binary squares and quadruple squares. They are illustrated in figures 2a and e and 3a and e. Notice that the scale of the axes is different in each case. Some physical properties of the models as well as the calculated solutions versus those extracted by K&H are shown in Table 1. It should be noted that owing to the dependence of the investigation depth on profile length in the K&H filtering method, the apparent current map shape is trapezoidal or triangular such that the maximum displayed depth for the pseudo-section is about 1/6 of the profile length (Ogilvy & Lee 1991). This is just because of the intrinsic behavior of the K&H filtering method and is not necessarily the depth of investigation. For instance, if the profile length is 6 km, the maximum recovered and displayed depth using K&H filtering is about 1 km, which is not reasonable for a VLF survey. The inversion method, on the other hand, does not have such a depth dependency on profile length and yields a rectangular section. In this case, the possibility of greater depth analysis, regardless of the profile length, is applicable. Although in the inversion method one can consider the depth of, say, 2 km for inversion model, it probably yields zero current density for that depth because approximately no anomalous body in the depth of 2 km can affect the VLF receiver at the surface. It should also be mentioned that the VLF-EM method itself is depth limited, theoretically being less than 50 m. However, in practical applications the range is slightly higher, up to 75–100 m, depending on the resistivity and thickness of the overburden. Figure 2. Open in new tabDownload slide (a) and (e) Schematic cross section of a dike and a synclinal structure, respectively; (b) and (f) in-phase component of H; (c) and (g) current density pseudo-section calculated using Karous–Hjelt filtering in conjunction with real solutions and spurious artifacts; (d) and (h) current density pseudo-section calculated from the inverse method with its real solutions. Figure 2. Open in new tabDownload slide (a) and (e) Schematic cross section of a dike and a synclinal structure, respectively; (b) and (f) in-phase component of H; (c) and (g) current density pseudo-section calculated using Karous–Hjelt filtering in conjunction with real solutions and spurious artifacts; (d) and (h) current density pseudo-section calculated from the inverse method with its real solutions. Figure 3. Open in new tabDownload slide (a) and (e) Schematic plan of a binary system of squares with the same resistivity but different dimensions and coordinates, and a quadruple system of squares with varying amounts of resistivity but fixed dimensions; (b) and (f) in-phase component of H for these two systems; (c) and (g) current density pseudo-section calculated from Karous–Hjelt filtering in conjunction with real solutions and spurious artifacts; (d) and (h) current density pseudo-section calculated from inverse method with its real solutions. Figure 3. Open in new tabDownload slide (a) and (e) Schematic plan of a binary system of squares with the same resistivity but different dimensions and coordinates, and a quadruple system of squares with varying amounts of resistivity but fixed dimensions; (b) and (f) in-phase component of H for these two systems; (c) and (g) current density pseudo-section calculated from Karous–Hjelt filtering in conjunction with real solutions and spurious artifacts; (d) and (h) current density pseudo-section calculated from inverse method with its real solutions. Table 1. Results from comparing the Karous–Hjelt filtering and inversion methods using synthetic models. The estimated locations from K&H or inversion methods are the locations of maximum current density. Location from K&H Model's cross section Horizontal extension [x1, x2] Vertical extension [z1, z2] Acceptable response [x, z] Spurious response [x, z] Location from inversion [x, z] Dipping Dike [400,480] [−10,−40] [435.3,−54.6] − [412.6,−35.2] Synclinal Structure [360,440] [10,35] [369.1,−29.9][400.7,−60][435.2,−29.9] − [362.7,−20.2][406.8,−64.9][437.3,−20] Binary Sq. (B) [60,70] [0,10] [64,−10] [103,−48] [66.7,−5] Sq. (A) [115,135] [10,30] [129.8,−25] [126.7,−19.9] Quadruple Sq. (1) [295,305] [10,20] [294.4,−20] [314,−35] [297.1,−14.8] Sq. (2) [335,345] [7.5,17.5] [339.6,−15.1] [339.6,−60] [338.4,−10.2] Sq. (3) [395,405] [35,45] Undetermined [374.2,−45] Undetermined Sq. (4) [495,505] [5,15] [499.1,−15.1] [455.5,−55.1] [497.2,−15.1] Location from K&H Model's cross section Horizontal extension [x1, x2] Vertical extension [z1, z2] Acceptable response [x, z] Spurious response [x, z] Location from inversion [x, z] Dipping Dike [400,480] [−10,−40] [435.3,−54.6] − [412.6,−35.2] Synclinal Structure [360,440] [10,35] [369.1,−29.9][400.7,−60][435.2,−29.9] − [362.7,−20.2][406.8,−64.9][437.3,−20] Binary Sq. (B) [60,70] [0,10] [64,−10] [103,−48] [66.7,−5] Sq. (A) [115,135] [10,30] [129.8,−25] [126.7,−19.9] Quadruple Sq. (1) [295,305] [10,20] [294.4,−20] [314,−35] [297.1,−14.8] Sq. (2) [335,345] [7.5,17.5] [339.6,−15.1] [339.6,−60] [338.4,−10.2] Sq. (3) [395,405] [35,45] Undetermined [374.2,−45] Undetermined Sq. (4) [495,505] [5,15] [499.1,−15.1] [455.5,−55.1] [497.2,−15.1] Open in new tab Table 1. Results from comparing the Karous–Hjelt filtering and inversion methods using synthetic models. The estimated locations from K&H or inversion methods are the locations of maximum current density. Location from K&H Model's cross section Horizontal extension [x1, x2] Vertical extension [z1, z2] Acceptable response [x, z] Spurious response [x, z] Location from inversion [x, z] Dipping Dike [400,480] [−10,−40] [435.3,−54.6] − [412.6,−35.2] Synclinal Structure [360,440] [10,35] [369.1,−29.9][400.7,−60][435.2,−29.9] − [362.7,−20.2][406.8,−64.9][437.3,−20] Binary Sq. (B) [60,70] [0,10] [64,−10] [103,−48] [66.7,−5] Sq. (A) [115,135] [10,30] [129.8,−25] [126.7,−19.9] Quadruple Sq. (1) [295,305] [10,20] [294.4,−20] [314,−35] [297.1,−14.8] Sq. (2) [335,345] [7.5,17.5] [339.6,−15.1] [339.6,−60] [338.4,−10.2] Sq. (3) [395,405] [35,45] Undetermined [374.2,−45] Undetermined Sq. (4) [495,505] [5,15] [499.1,−15.1] [455.5,−55.1] [497.2,−15.1] Location from K&H Model's cross section Horizontal extension [x1, x2] Vertical extension [z1, z2] Acceptable response [x, z] Spurious response [x, z] Location from inversion [x, z] Dipping Dike [400,480] [−10,−40] [435.3,−54.6] − [412.6,−35.2] Synclinal Structure [360,440] [10,35] [369.1,−29.9][400.7,−60][435.2,−29.9] − [362.7,−20.2][406.8,−64.9][437.3,−20] Binary Sq. (B) [60,70] [0,10] [64,−10] [103,−48] [66.7,−5] Sq. (A) [115,135] [10,30] [129.8,−25] [126.7,−19.9] Quadruple Sq. (1) [295,305] [10,20] [294.4,−20] [314,−35] [297.1,−14.8] Sq. (2) [335,345] [7.5,17.5] [339.6,−15.1] [339.6,−60] [338.4,−10.2] Sq. (3) [395,405] [35,45] Undetermined [374.2,−45] Undetermined Sq. (4) [495,505] [5,15] [499.1,−15.1] [455.5,−55.1] [497.2,−15.1] Open in new tab Next, each model is explained individually. 3.1. Dipping dike A dipping dike with 10 m thickness, spanning horizontally from 400 m to 480 m and vertically from −10 m to −40 m is assumed as a model (figure 2a). The in-phase component of VLF data due to this model is plotted in figure 2b. The current density contour map relating to the dipping dike, recovered using the K&H approach, is open ended (figure 2c). The maximum point of the contour map is detected at (435.3, −54.6). The recovered current density contour map using the inversion method is also displayed in figure 2d and the current contours in this condition turn out to be closed at a shallower depth. Although the contour map is not indicative of the true dip direction of the dike for either of the methodologies (Ogilvy & Lee 1991), asymmetry gives an overall rewarding clue to the general dip direction. The maximum current point resulting from inverse method lies in (412.6, −35.2). This point is closer to the upper surface of the dike where the highest amount of current, in reality, is induced (figure 2d). Accordingly, when compared with K&H filtering, the inversion method approximates the real current behavior more precisely. 3.2. Synclinal structure A synclinal structure, extending horizontally from 400 to 480 m and vertically from −10 to −40 m is assumed as our model (figure 2e). The real component of the VLF data corresponding to this model is shown in figure 2f. Considering this structure, the accuracy and superiority of our method relative to that of K&H, specifically for shallower bodies, is well revealed. Implementation of the K&H filter to this structure demonstrates three high current spots (figure 2g and Table 1). This can definitely be confusing for an interpreter as to whether these points show three separate masses with the most conductive one laying in the middle, or a deeper location, a synclinal structure or two shallow bodies with a spurious artifact emerging between them in a lower position. Whatever the interpretation could be, these maxima do not illustrate the true vertical location of the synclinal structure (figure 2g and Table 1). The upper surface of both arms and their joint locations are estimated to be far deeper than their real coordinates. In contrast, the inversion approach pinpoints the upper side of the arms relatively precisely and marks a steady and uniform high current zone corresponding to the apparent shape of the syncline (figure 2h and Table 1). This area is deeper than the syncline body, however. Just like the K&H filter, the inverse method does not yield a realistic sign of the lower part of the anomalous structure, although at least it shows that the vertical extent is limited. Even in horizontal aspects, the inverse methodology acts very well, and the lateral edges of the arms are delineated acceptably. Nevertheless, these points in K&H methodology lie somewhere toward the inner part of the arms and lower than their true locations (figure 2g and h). 3.3. Binary squares The purpose of examining two neighboring anomalous bodies with the same resistivity and different dimensions and coordinates is to analyze the effect of interference. Figure 3a shows two bodies with rectangular sections as well as their known horizontal and vertical extensions. Figure 3b shows the in-phase VLF data for models shown in figure 3a. Two maxima were obtained by the K&H method at (129.8, −25) and (64, −10) for squares A and B, respectively (figure 3c). These spots are within the bodies, and could be helpful in an overall estimation of the mass locations, but they are not exact. The real high current positions should occur at the uppermost side of the squares, whereas the maximum current points for squares A and B, obtained by inversion method, are more exact so that the peaks are at (126.7, −19.9) and (66.7, −5), in order (figure 3d). Another point of critical importance is the appearance of spurious artifacts in the K&H approach when dealing with interfering sources. As shown in figure 3c, a high current zone appears at (103, −48) that is completely misleading. This could be deceptive and lead to incorrect deduction from the survey unless further geophysical and geological information is provided to discriminate the real effects from unreal ones. An advantage of our methodology over K&H filtering is elimination of spurious high current zones. Therefore, while providing the current maps using the inversion method, there is a certainty that no current maximum point is fallacious. 3.4. Quadruple squares In this section, four squares 1, 2, 3 and 4 of the same size but different resistivities and depths are assessed. The locations and extensions of these models are shown in figure 3e, while calculated VLF data from them is plotted in figure 3f. The resistivity of square 2 is 10 times greater than the resistivity of squares 1 and 3, and also 20 times greater than the resistivity of square 4. So, its current effect might be dominated by its neighbors, in particular square 1. Applying K&H filtering, three high spots are disclosed (figure 3g and Table 1). It could be seen that even though the lower surface of the targets is detected, the high current zone of square 3 is not observable. This may be found in a far deeper location if further depth examination is processed via tapering the data (lengthening the profile) (Ogilvy & Lee 1991). Nonetheless, if this is done, because of numerous spurious artifact appearances and a lack of ability to recognize square 3 effects among all those unreal high current points, this initiative is absurd in practice. On the contrary, the inversion method could give an acceptable estimate of the targets (figure 3d and Table 1), and show their center or upper sides. As is obvious, in this condition, there is no spurious effect to confuse the interpreter. However, the current effect of square 4 is not yet observed. This is definitely because of the field domination of its nearby masses. 4. Real data ‘Ghoori Ghale’ cave is 3140 m long; the longest cave in Iran. It is located about 87 km northwest of Kermanshah city in the west of Iran. Its entrance coordinate is 34˚53.98′ N, 46˚30.11′ E. According to Iranian tectonics, the cave belongs to the Zagros thrust zone and the main trends in the area are completely related to north Zagros zone and its evolution. The genesis of the area goes back to the Mesozoic era, so most of its rocks are from the Jurassic and Cretaceous. Surface rocks, having numerous outcrops, contain green-red tuffaceous siltstone and sandstone, red banded chert, limestone, mudstone and radiolarians (figure 4). Tension from the Arabian plate to Iranian one is the main reason for the faulting of the latter plate and as a result the Zagros Mountains have been formed. Most of the existing faults in the area are in the parallel direction to the Zagros Mountains, running from northwest to southeast. Faulting in rocks and then transmission of water through them have caused erosion of sandstones and limestones, so the cave has gradually been created. Besides the main cave, some galleries that are not very long have been created, more and less, in the vicinity of the cave. Figure 4. Open in new tabDownload slide Geological map (Sadeghian & Delavar 2006) containing the study area (rectangle in the right side), as well as a Google map of the Ghoori Ghale cave and data profiles on the left side. Figure 4. Open in new tabDownload slide Geological map (Sadeghian & Delavar 2006) containing the study area (rectangle in the right side), as well as a Google map of the Ghoori Ghale cave and data profiles on the left side. The field work was done using GSM-19 apparatus at a frequency of 19.6 kHz. Three profiles were conducted at the strike approximately perpendicular to the likely cave route (figure 4). Here, we have chosen the data from profile No. 2 with a length of 210 m. Real and imaginary components of VLF data, acquired over profile No. 2, as well as current density pseudo-sections obtained using K&H filtering and inversion method are, respectively, displayed in figure 5. The local minimums of current density, which are representatives of high resistivity bodies, are indicated with red circles. Both the K&H and inversion method have discriminated the presence of high resistivity zones but there are just three minimum current density circles in the inverted map against seven circles in the K&H map. So, discrimination of global minimums is more probable in the inversion method than the K&H one. Besides, the inversion method yields a full current density map of the subsurface beneath the profile line while the K&H method fails to recover the lower left and right sides of the map. Therefore, according to the goal of this survey, the result was that there is a big resistive zone, which in our case is a cave, at a distance of 100 m from beginning of the profile that at its center has a depth of ≈32 m. Furthermore, a small canal or cave is observable at a distance of 205 m from the profile origin, which at its center is ≈16 m in depth. Figure 5. Open in new tabDownload slide (a) Real and imaginary components of H relating to profile 1 over the Ghoori Ghale cave; (b) current density pseudo-section calculated from Karous–Hjelt filtering along with all the calculated solutions that, in terms of being realistic or unrealistic, are not distinguishable from each other and (c) current density pseudo-section calculated from the inverse method with its solutions. Figure 5. Open in new tabDownload slide (a) Real and imaginary components of H relating to profile 1 over the Ghoori Ghale cave; (b) current density pseudo-section calculated from Karous–Hjelt filtering along with all the calculated solutions that, in terms of being realistic or unrealistic, are not distinguishable from each other and (c) current density pseudo-section calculated from the inverse method with its solutions. 5. Results and discussion After comparing the two mentioned methods, some characteristics, upsides and downsides of both procedures can be recognized as follows: In both methods, the depth estimation is erratic specifically when the body is deep. In the case of shallow anomalies, nevertheless, the precision of the inverse method is far more trustworthy and surpasses the one offered by K&H. Even in the horizontal aspect, the inverse method excels beyond the K&H method and pinpoints the lateral spread of the underground masses more realistically. When using the K&H filtering method for interfering sources, due to spurious artifacts, making any statement about deeper high current zones is undetermined and accompanied by doubt and uncertainty, while this issue is removed in inverse methodology and, due to the nature of this procedure, no artifact appears. Investigating higher depths in the K&H approach, the width of the pseudo-current map decreases and a trapezoidal or triangular map is formed. Therefore, some precious and vital data might be missing, whereas in our inverse-based theory, a full rectangular pseudo-current map is attained. Concerning uniformly conductive synclinal structures, current maps pertaining to the inverse method track the apparent shape of an underground anomalous body, unlike that of K&H filtering, only showing three high current spots, which raise doubt (figure 2g and h). Regarding the vertically elongated bodies like dikes, neither of the approaches could determine the true dip angle, yet they give a whole insight into the dip direction. Even though the K&H method may have a more powerful distinguishing property than the inverse method, when dealing with a multitude of separate neighboring masses at different depths, this characteristic is not practically functional because of the advent of spurious artifacts. 6. Conclusion In this study, the current density pseudo-section, related to VLF data over a profile, is mapped using a damped least squares method. Even though the new methodology still miscalculates the vertical coordinate of deep anomalies, in the case of shallow targets it works very well and is superior to the Karous–Hjelt filtering method both in vertical and horizontal location estimations. Moreover, for multiple close targets, the appearance of spurious artifacts that is the main drawback of Karous–Hjelt filtering is not observed in the inverse theory. Accordingly, vagueness and uncertainty in the interpretation process is removed to a great extent. For interpretation of real data, the inversion method can recover the global minimum (or maximum) with more accuracy than the well-known K&H filtering method. Here, the horizontal and vertical location of the Ghoori Ghale cave is better recovered than by using the K&H method. Acknowledgements Conflict of interest statement. None declared. References Al Bulushi A.M. , Al Wardi M. , Al Shaqsi B. , Sundararajan N. , 2016 . Mapping of subsurface fault structures by VLF-EM method in Al Khoud area, Muscat, Sultanate of Oman , Arabian Journal of Geosciences , 9 , 349 . Google Scholar Crossref Search ADS WorldCat Aster R.C. , Borchers B. , Thurber C. , 2003 . Parameter Estimation and Inverse Problems , Elsevier Publications . Google Preview WorldCat COPAC Fraser D.C. , 1969 . Contouring of VLF-EM data , Geophysics , 34 , 958 – 967 . Google Scholar Crossref Search ADS WorldCat Gharibi M. , Pedersen L.B. , 1999 . Transformation of VLF data into apparent resistivities and phases , Geophysics , 64 , 1393 – 1402 . Google Scholar Crossref Search ADS WorldCat Gnaneshwar P. , Shivaji A. , Srinivas Y. , Jettaiah P. , Sundararajan N. , 2011 . Very-low-frequency electromagnetic (VLF-EM) measurements in the Schirmacheroasen area, East Antarctica , Polar Science , 5 , 11 – 19 . Google Scholar Crossref Search ADS WorldCat Karous M. , Hjelt S.E. , 1977 . Determination of apparent current density from VLF measurements . Report, Department of Geophysics, University of Oulu, Finland. Contribution No. 89 , 1 – 18 . Google Preview WorldCat COPAC Karous M. , Hjelt S.E. , 1983 . Linear filtering of VLF dip-angle measurements , Geophysics Prospect ., 31 , 782 – 794 . Google Scholar Crossref Search ADS WorldCat Khalil M.A. , Abbas M.A. , Monterio Santos F.A. , Mesbah H.S.A. , Massoud U. , 2010 . VLF-EM study for archaeological investigation of the labyrinth mortuary temple complex at Hawara area, Egypt , Near Surface Geophysics , 8 , 203 – 212 . Google Scholar Crossref Search ADS WorldCat Khalil M.A. , Monterio Santos F.A. , 2010 . Comparative study between filtering and inversion of VLF-EM profile data , Arabian Journal of Geosciences , 4 , 309 – 317 . Google Scholar Crossref Search ADS WorldCat McNeill J.D. , Labson V.F. , 1991 . Geological mapping using VLF radio fields , In: Electromagnetic Methods in Applied Geophysics, 2 (applications), part B, Investigations in Geophysics, Society of Exploration Geophysicists , pp. 521 – 640 , ed. Nabighian M.N . Tulsa, Oklahoma . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Monteiro Santos F.A. , Mateus A. , Figueiras J. , Gonçalves M.A. , 2006 . Mapping groundwater contamination around a landfill facility using the VLF-EM method – a case study , Journal of Applied Geophysics , 60 , 115 – 125 . Google Scholar Crossref Search ADS WorldCat Ogilvy R.D. , Lee A.C. , 1991 . Interpretation of VLF-EM in-phase data using current density pseudosections , Geophysics Prospect ., 9 , 567 – 580 . Google Scholar Crossref Search ADS WorldCat Parker M.E. , 1980 . VLF electromagnetic mapping for strata-bound mineralization near Aberfeldy, Scondland , Trans. Inst. Min. Metall., Sect B , 89 , B123 – B133 . WorldCat Paterson N.R. , Ronkav V. , 1971 . Five years of surveying with the very low frequency electromagnetic method , Geoexploration , 9 , 7 – 26 . Google Scholar Crossref Search ADS WorldCat Philips W.J. , Richards W.E. , 1975 . A study of the effectiveness of the VLF method for the location of narrow mineralized fault zones , Geoexploration , 13 , 215 – 226 . Google Scholar Crossref Search ADS WorldCat Ramesh Babu V. , Ram S. , Srinivas R. , Veera Bhaskar D. , Bhattacharya A.K. , 2004 . VLF-EM surveys for uranium exploration in Dulapali area, Raigarh district, Madhyapradesh, India , Journal of Geophysics , 5 , 27 – 33 . WorldCat Ramesh Babu V. , Ram S. , Sundararajan N. , 2007 . Modeling of magnetic and VLF-EM with an application to basement fracturesda case study from Raigad India , Geophysics , 71 , 133 – 140 . Google Scholar Crossref Search ADS WorldCat Reitz J.R. , Milford F.J. , 1966 . Foundations of Electromagnetic Theory , Addison Wesley Publications , Tokyo . Google Preview WorldCat COPAC Sadeghian M. , Delavar S.T. , 2006 . 1:100000 Geological Map of Kamyaran , Geological Survey of Iran publications . Google Preview WorldCat COPAC Sedelnikov E.S. , Gridchin V.M. , Gardev S.G. , Lunin Yu G. , Rogachev B.V. , 1971 . Method of Long Wavelength Radiokip with SDVR-3 Equipment (in Russian) , Tchinigri Moscow . Google Preview WorldCat COPAC Sinha A.K. , Hayles J.C. , 1988 . Experiences of local loop VLF transmitter for geological studies in the Canadian Nuclear Fuel Waste Management program , Geoexploration , 25 , 37 – 60 . Google Scholar Crossref Search ADS WorldCat Sundararajan N. , Narasimah Chary M. , Nandakumar G. , Srinivas Y. , 2007 . VES and VLF-an application to groundwater exploration Khammam India , The Leading Edge , 26 , 708 – 716 . Google Scholar Crossref Search ADS WorldCat Sundararajan N. , Ramesh Babu V. , Chatruvedi A.K. , 2011 . Detection of basement fractures favorable to uranium mineralization from VLF-EM signals , Journal of Geophysics and Engineering , 8 , 330 – 340 . Google Scholar Crossref Search ADS WorldCat Sundararajan N. , Ramesh Babu V. , Shiva Prasad N. , Srinivas Y. , 2006 . VLFPROS – A MATLAB code for processing of VLF-EM data , Computers & Geosciences , 32 , 1806 – 1813 . Google Scholar Crossref Search ADS WorldCat Tikhonov A.N. , Arsenin V.Y. , 1977 . Solution of Ill-Posed Problems , Winston & Sons Washington . Google Preview WorldCat COPAC Timur E. , 2012 . VLF-R studies in the Agora of Magnesia archaeological site, Aydin, Turkey , Journal of Geophysics and Engineering , 9 , 697 – 671 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2020. Published by Oxford University Press on behalf of the Sinopec Geophysical Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - A one-step linear inversion to recover current density pseudo-section from VLF-EM signals
JF - Journal of Geophysics and Engineering
DO - 10.1093/jge/gxz063
DA - 2020-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/a-one-step-linear-inversion-to-recover-current-density-pseudo-section-TP39Jhp9LP
SP - 1
VL - Advance Article
IS - 
DP - DeepDyve
ER -