\documentclass[twoside]{article}
\usepackage[silent,pdf,hyper]{elxshort}
\usepackage{color}
\usepackage{hyperref}
\journalid{ES}{19}{3}{19}
\citnumber{}{3}
\paperid{2019ES000664}
\papertype{RAR}
\cpright{}{2019}
\lefthead{aleshin and malygin}
\righthead{approach to inter-well radio wave}
\received{13}{April}{2019}
\accepted{16}{May}{2019}
\published{22}{May}{2019}

\authaddr{I. M. Aleshin, Schmidt Institute of Physics of the Earth RAS, Moscow, Russia (ima{@}ifz.ru)}

\begin{document}

\title{Machine learning approach to inter-well radio wave survey data imaging}

\author{1}{Aleshin}{I. M.}{}{1}{2}{}
\author{1}{Malygin}{I. V.}{}{1}{}{}

\affil{1}{Schmidt Institute of Physics of the Earth RAS}{Moscow}{}{Russia}
\affil{2}{Laverov Federal Center for Integrated Arctic Research RAS}{Arkhangelsk}{}{Russia}

\keywd{inter-well scanning}
\keywd{radio wave survey}
\keywd{machine learning}
\keywd{kNN-algorithm}


\abstract{Inter-well measurements are used to reduce drilling costs with no reduce small kimberlite body detection. The radio wave method enables measurement of the apparent absorption coefficient that is proportional to the effective electrical resistance of the rock. Our point is to build a three-dimensional model of distribution of electrical properties of inter-well space throughout the entire exploration region.  The measured data is distributed unevenly because data points are grouped along the linear clusters. The distance between neighbor points composing a cluster is much smaller than distance between clusters. In terms of geostatistics, this means a significant spatial anisotropy of data distribution that is difficult to take into account using standard geostatistical approach. We have shown that the problem could be solved by methods developed within the theory of machine learning. To build a three-dimensional model of attenuation coefficient we used a modified method of $k$-nearest neighbors.}

\section{Introduction}

Currently, deposits of diamonds directly accessible from the surface are almost exhausted in the Western Yakutia. Search for kimberlite bodies is carried out in areas where traditional geological and geophysical studies are ineffective [\itc{Shmakov,} \reflink{shmak17}{2017}]. For areas covered by sedimentary rocks, as well as traps, the only direct method to search for kimberlites is to run a network of borehole measurements. To reduce the cost of work it is desired to increase the distance between wells. However, this increases the risk of missing small kimberlite pipes. To avoid this, cross-well survey methods are used, in particular, radio wave methods.
%%[\itc{Kevorkyanc at al.,}\reflink{kev05}{2005}].
The radio wave scanning technique was developed in the middle of the last century (see, for review, [\itc{Petrovskij,} \reflink{pet71}{1971}]) and is being actively used today [\itc{Istratov at al.,} \reflink{ist06}{2006}]. This technique is used not only when searching for kimberlite pipes [\itc{Tolstov at al.,} \reflink{tol18}{2018}], ore [\itc{Kuznetsov,} \reflink{kuz08}{2008}] and oil [\itc{Istratov at al.,} \reflink{ist00}{2000}] deposits, etc., but also for natural hazards [\itc{Cherepanov,} \reflink{cher17}{2017}] and technological processes monitoring [\itc{Istratov at al.,} \reflink{ist09}{2009}].\vspace*{2ex}

The idea of the method is to estimate the attenuation of an electromagnetic wave as it passes between two wells. The source and receiver of the electromagnetic field is placed in adjacent wells to measure the attenuation of the electric field amplitude. Lower rock resistance corresponds to higher radio wave absorption. Therefore, the absorption coefficient at a fixed frequency is proportional to the electrical conductivity of the medium [\itc{Petrovskij,} \reflink{pet71}{1971}]. There are two kinds of inter-well scanning techniques. First one is named ``fan'' method shown on panel A of \figref{1}.  The position of the emitter in the well varies with a given step. At each position, the receiver fixes the amplitude of the electric field emitted by the source. After that, the source moves to the next position. The source also moves throughout the entire working interval of the well. Such a measurement scheme allows to obtain a detailed picture of the electrical properties of the inter-well space in the plane passing through both wells. The joint interpretation of the set of sections obtained in this way allows to obtain a three-dimensional image of the electrical properties of the medium [\itc{Kuznetsov,} \reflink{kuz12}{2012}].

\begin{figure}[t]          % Fig 1
\figurewidth{20pc}
\setimage{}{1.3ex}{20pc}{}{2019es000664-f01}
\caption{Measurement schemes used in radio wave research. A -- fan schema, the position of the source is fixed, the receiver moves through the entire working range with a given step. Then the source is shifted to the next position and the measurements are repeated, etc. B -- synchronous measurements, the source and receiver are moved along the well at the same time.}
\end{figure}

If fan survey was chosen $n^{n}$ measurements must be made where $n$ is the number of stops (measurement positions) per well. In practice the survey is often limited to synchronous dipping of the source and receiver in the neighboring wells. In this case, the number of measurements is equal to $n$, which significantly reduces the quantity of measurements, but at the same time the information content of performed measurements also decreases. In fact, such  measurement scheme allows us to get only the average value of the apparent attenuation coefficient corresponding to the midpoint of the line connecting the source and receiver positions. Thus the tomographic method becomes inapplicable. To build a three-dimensional model of the medium, we need to use an alternative interpolation procedure. [\itc{Aleshin and Zhandalinov,} \reflink{al09}{2009}]  proposed to construct horizontal sections of the model using kriging [\itc{Isaaks and Srivastava,}\reflink{is89}{1989}]. These parameters should lead to an interpolation error, which is comparable to the accuracy of the input data. However, while constructing a real three-dimensional model, the direct application of the regression methods is impossible because of the extreme anisotropy of the distribution of input data. In addition, even in the two-dimensional case, the solution turns out to be too smooth while our point is to increase image contrast.

One possible alternative is to use machine learning methods. In this case, the interpolation procedure is based on a specific analysis of the source data, which is called training. In machine learning, our task can be classified as the analysis of ordered data. The traditional approach is to build models based on deep neural network architectures, for example, convolutional or recurrent [\itc{Nikolenko at al.,} \reflink{nik18}{2018}]. However, the data we have is not enough for quality training of such models. We used the simplest implementation of the $k$-Nearest Neighbors algorithm (abbreviated kNN). This is known as simple but effective method of data analysis. The kNN method belongs to the class of so-called lazy algorithms. Such algorithms do not require a long preliminary training. The decision rules are based on calculating the distances between the data objects [\itc{Zhuravlev et al.,} \reflink{jur06}{2006}]. Actually, training is reduced to the calculation of the distance matrix from a given point to all input data.

\section{Methods}

In this work, we use the data of ALROSA company (\href{http://www.alrosa.ru}{http://www.alrosa.ru}). Synchronous inter-well radio wave measurement  was performed at one of sites in Yakutia. We are unable to use additional data related to the site. Therefore, we restricted our effort to the construction of a three-dimensional model of conductive properties. In later parts, we will discuss only model of spatial distribution of the apparent absorption coefficient of the medium between the wells. To analyze the inter-well measurement data, it is necessary to select a model of the propagation of the electric field emitted by the source. Since our data is obtained from synchronous measurements we can estimate only the average value of the absorption coefficient of the medium along the straight line connecting the source and the receiver.  Therefore, we confine ourselves to the simplest model of wave propagation, described by the formula of the radiation field of an electric dipole in a homogeneous isotropic medium

$$
E = E_0 \exp(-q/R) / R \sin\theta.
$$

\noindent Here, $E$ is the polar component of the electric field, $E_{0}$ is the amplitude of the emitted wave, $R$ is the distance between the source and receiver positions. With synchronous measurements, the source and receiver are placed approximately at the same depth, so the polar angle can be set equal to $\pi/2$. Then the absorption coefficient is

$$
q = - \ln \left(RE/E_0\right).
$$

Let $Q=\{q_n\}$ denote the set of input data: the values of the apparent damping factor measured at $N$ points with coordinates $\vec r_n=\{x_{n}, y_{n}, z_{n}\}$. Here, the $x, y$ coordinates determine the position of a point in the horizontal plane, $z$ is the depth measured from sea level. In the kNN algorithm, the value of $q$ at an arbitrary point $\vec r = \{x,y,z\}$ should be  calculated by the formula

\begin{equation}\label{1}
q\left(\vec r\right)=\sum_{k=1}^K{w_k\left( \vec r, \vec r_k\right) q_k}, \quad
\sum_{k=1}^K w_k = 1,
\end{equation}

\noindent the summation is over $K$ point closest to $\vec r$.  The number $K$ is a free parameter of the algorithm (hyperparameter), which requires additional determination. Instead of the number of neighbors, the radius of the sphere with center at the point $\vec r$ is used as the hyperparameter and all points inside the sphere are considered as neighbors. It is natural to choose the Euclidean metric as the distance between points:

\begin{equation}\label{2}
R\left(\vec r_1, \vec r_2 \right) =
\sqrt{\left(x_1 - x_2\right)^2+\left(y_1 - y_2\right)^2+\left(z_1 - z_2\right)^2}.
\end{equation}

\noindent The weight function $w_k\left(\vec r, \vec r_k\right)$ depends on the distance between the current point $\vec r$ to the corresponding point with a given value. As a weight function, a value inversely proportional to the distance is usually used:
$$
w_k \left(\vec r, \vec r_k\right) \sim 1/R\left(\vec r, \vec r_k \right).
$$
Without taking into account the spatial location of the points of weight are the same for all members of $w_k=1/K$.

\begin{figure}         % Fig 2
\figurewidth{20pc}
\setimage{}{}{20pc}{}{2019es000664-f02}
\caption{The effect of the scale factor $\lambda$ on the nearest point selection.  Anisotropy of data distribution leads to nearest neighbors are within the same measurement group for any point (A). Scaling the horizontal axes corrects this situation (B).}
\end{figure}


The distribution of data obtained by the radio wave scanning method is strongly anisotropic. The depth step is 5 m with a well length of the order of 500 m while the distance between the nearest wells is approximately 200 m (see \figref{2}). We cannot use classical methods of geostatistics [\itc{Isaaks and Srivastava,} \reflink{is89}{1989}] to build a three-dimensional model of the media. Moreover, the application of the method of nearest neighbors also requires its modification. To reduce the difference in horizontal and vertical scales, we will redefine the metric \eqref{2} by entering a dimensionless scale factor $\lambda$:
\begin{multline}\label{3}
R\left(\vec r_1, \vec r_2 \right) = \\
\sqrt{
   \left(x_1 - x_2\right)^2 / \lambda^2+
   \left(y_1 - y_2\right)^2  /\lambda^2 +
   \left(z_1 - z_2\right)^2
}.
\end{multline}

One can expect that a proper choice of the parameter value would compensate the anisotropy of the data, but we have no criterion for making this choice. Therefore, the scale factor, along with the number of neighbors $K$, is another hyperparameter of the problem. We used the cross-validation to determine the hyperparameters. This approach, along with the hold-out method, is standard in machine learning theory. The source data is divided into $M$ groups ($M=5$), each of these groups is used for testing. In our case, to estimate the quality of the solution, the total rms value of the difference between the deferred values of the observed data, calculated from the remaining data using the formula \eqref{1} with the metric \eqref{3}, is used:

$$
\kappa(K,\lambda) = 1/M \sum_{m=1}^M{\kappa^{(m)}(K,\lambda)}
$$

$$
\kappa^{(m)} (K,\lambda) =
$$

$$
1 - \sum_{i=1}^{N/M} \left( q_i^{(m)} - q(\vec r_i; K,\lambda) \right)^2 / \sum_{i=1}^{N/M} \left(q_i^{(m)} - \bar q^{(m)}\right)^2
$$

$$
 \bar q^{(m)} = \sum_{i=1}^{N/M} q_i^{(m)}.
$$

\noindent Interpolant $q(\vec r_i; K,\lambda)$ is calculated by the formula \eqref{1} using the metric \eqref{3} but the holdout  data are not taken into account. A coefficient of determination was calculated on a 25 by 25 grid  with a unit step for both parameters. The result of the calculations is shown in \figref{3}. The coefficient distribution has the shape that is typical for multi-parameter optimization problem. To select the hyperparameter values we have fixed the level 0.7 on the determination coefficient map. It approximately  corresponds to the 80-percent correlation between the model and the input data. Intersection of the median of the triangle, which is formed by the coordinate axes and the straight line approximation of the 0.70 level of the determination coefficient to select of the hyperparameter values equal $K=11$ and $\lambda=10$~m.

\begin{figure}        % Fig 3
\figurewidth{20pc}
\setimage{}{}{20pc}{}{2019es000664-f03}
\caption{The determination coefficient distribution calculated on the grid of parameters.  An acceptable values of  the hyperparameters corresponds to area where determination coefficient value exceeds level 0.7. To build the model, the values determined by the intersection point of the median and hypotenuse of the triangle formed by the axes of coordinates and the straight line approximating the level 0.7 are chosen.}
\end{figure}


\begin{figure*}           % Fig 4
\figurewidth{35pc}
\setimage{}{}{35pc}{}{2019es000664-f04}
\caption{Horizontal (panels B, D) and vertical (panel C) cross-sections of the model. Depth axis origin corresponds to the sea level, the position of the horizontal axes is consistent with the geometry of the site. The location of the vertical section corresponds to the line $Y=0$ (line on panel A). Horizontal cross-sections are shown on panels B, D correspond to depths $Z=-560$ and $Z=-250$ meters (black and white dotted lines on the panel C, correspondingly ).}
\end{figure*}

Since the values of the hyperparameters are determined we can construct
the attenuation coefficient image. The calculation of horizontal and vertical cross-sections is implemented in the Python programming language (\href{https://www.python.org}{https://www.\linebreak python.org}) using the Scikit-learn package collection (\href{https://scikit-learn.org}{https://scikit-learn.org}).

\figref{4} shows one vertical and two horizontal cross-sections of the model. It can be seen that the model constructed allows to localize objects whose horizontal size are significantly smaller than the distance between the wells. As an example, the areas of high attenuation coefficient values located at a depth of -560 m and horizontal coordinates $X=2950$, $X=4750$ and $X = 5150$ meters can be cited. Their positions are indicated by the corresponding dashed lines on panels A and C, \figref{4}.


\begin{figure*}      % Fig 5
\figurewidth{35pc}
\setimage{}{}{35pc}{}{2019es000664-f05}
\caption{Vertical cross-sections of the model correspond to the planes $X=2950, X=4750$ and $X=5150$ meters. In the previous figure, their positions are marked with vertical dotted lines. In the figures, compact areas of increased values of attenuation coefficient are clearly visible. Horizontal size of the arias  are about of 100 m.}
\end{figure*}

For clarity, the \figref{5} shows the vertical cross-sections corresponding to these lines, on which the corresponding areas are also clearly visible.

\section{Conclusions}

Inter-well synchronous survey data can be used to construct a three-dimensional model of inter-well medium conductivity using the kNN method. The strong anisotropy of the input data can be reduced by modifying the spatial metric, which determines the distance between the data. This can be achieved by introducing a scaling factor that scales in the horizontal direction. The approach used makes it possible to obtain a fairly contrasting image of inhomogeneous areas. In particular, it gives the possibility to make the outlining of areas with linear size smaller than the distance between the wells. The model building process does not depend on the physical model used to interpret the measurement. Refining the model of the physics of wave propagation between wells will improve the quality of image construction without the change in image construction procedure. Of course, the model can be improved by drawing additional data (geological, seismic, magnetic) for their joint interpretation.

\acknowledgments{This work was supported by budget projects of the Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, and the Laverov Federal Center for Integrated Arctic Research, Russian Academy of Sciences (project number AAAA-A18-118012490072-7).}

\begin{references}

\reference{al09}
\reftype{2}
\refauthf{Aleshin}{I. M.}{}
\refauth{Zhandalinov}{V. M.}{}
\refyear{2009}
\reftitle{Application of interpolation procedures for presentation of data electromagnetic wave lightning}
\refjour{Russ. J. Earth Sci.}
\refvol{11}
\refnumb{1}
\refpages{1--4}
\refdoi{10.2205/2009ES000430}
\endreference{}

\reference{cher17}
\reftype{1}
\refauthf{Cherepanov}{A. O.}{}
\refyear{2017}
\reftitle{Multi-frequency radio wave measurements in wells to monitor
the process of thawing MMP (example of the Russkoe oil field, Western
Siberia)}
\refjour{KRAUNZ Bulletin. Series: Earth Science}
\refvol{4}
\refpages{118--123}
\endreference{in Russian}

\reference{is89}
\reftype{3}
\refauthf{Isaaks}{E. H.}{}
\refauth{Srivastava}{R. M.}{}
\refyear{1989}
\refbook{Applied Geostatistics}
\refpages{589}
\refpubl{Oxford University Press}
\reflocat{New York}
\endreference{}

\enlargethispage{-14\baselineskip}

\reference{ist00}
\reftype{1}
\refauthf{Istratov}{V. A.}{}
\refauth{Lysov}{M. G.}{}
\refauth{Chibrikin}{I. V.}{}
\refgroup{et al.}
\refyear{2000}
\reftitle{Radio wave geointroscopy (RWGI) of inter-well space in oil fields}
\refjour{Geophysics}
\refvol{Special issue}
\refpages{59--68}
\endreference{in Russian}

\reference{ist06}
\reftype{1}
\refauthf{Istratov}{V. A.}{}
\refauth{Skrinnik}{A. V.}{}
\refauth{Perekalin}{S. O.}{}
\refyear{2006}
\reftitle{New equipment for radio wave geointoscopy of rocks in the
interwell space ``RWGI-2005''}
\refjour{Instruments and systems for exploration geophysics}
\refnumb{1}
\refpages{37--43}
\endreference{in Russian}

\reference{ist09}
\reftype{1}
\refauthf{Istratov}{V. A.}{}
\refauth{Kolbenkov}{A. V.}{}
\refauth{Perekalin}{E. V.}{}
\refauth{Lyax}{S. O.}{}
\refyear{2009}
\reftitle{Radio \ wave \ monitoring \ method \ of \ technological \ processes \ in \ the \ interwell \ space}
\refjour{KRAUNZ Bulletin. Series: Earth Science}
\refvol{14}
\refpages{59--68}
\endreference{in Russian}

\reference{kev05}
\reftype{1}
\refauthf{Kevorkyanc}{S. S.}{}
\refauth{Abramov}{V. Y.}{}
\refauth{Kovalev}{Y. D.}{}
\refyear{2005}
\reftitle{Well radio wave complex for searching for kimberlite pipes in
Western Yakutia}
\refjour{Geophysics}
\refvol{3}
\refpages{56--64}
\endreference{in Russian}

\reference{kuz08}
\reftype{1}
\refauthf{Kuznetsov}{N. M.}{}
\refyear{2008}
\reftitle{Experience of the radiowave geointroscopy of the interwell
space for the exploration of a gold-copper deposit}
\refjour{Exploration and protection of mineral resources}
\refnumb{12}
\refpages{27--29}
\endreference{in Russian}

\reference{kuz12}
\reftype{1}
\refauthf{Kuznetsov}{N. M.}{}
\refyear{2012}
\reftitle{The 3D method of processing the data of radio wave scanning of
the interwell space}
\refjour{KRAUNZ Bulletin. Series: Earth Science}
\refnumb{1}
\refpages{240--246}
\endreference{in Russian}

\reference{nik18}
\reftype{3}
\refauthf{Nikolenko}{S. I.}{}
\refauth{Kadurin}{A. A.}{}
\refauth{Arxangel'skaya}{E. O.}{}
\refyear{2018}
\refbook{Deep learning}
\refpages{479}
\refpubl{Piter}
\reflocat{St. Petersburg}
\endreference{in Russian}

\reference{pet71}
\reftype{3}
\refauthf{Petrovskij}{A. D.}{}
\refyear{1971}
\refbook{Radio wave methods in underground geophysics}
\refpages{224}
\refpubl{Nedra}
\reflocat{Moscow}
\endreference{in Russian}

\reference{tol18}
\reftype{4}
\refauthf{Tolstov}{A. V.}{}
\refauth{Zinchuk}{N. N.}{}
\refauth{Serov}{I. V.}{}
\refyear{2018}
\reftitle{Main results of research and experimental-methodical works of
``ALROSA'' (PJSC)}
\refbook{Efficiency of geological exploration for diamonds: forecasting
and resource, methodical, innovative and technological ways to increase it}
\refpubl{``ALROSA'' company}
\reflocat{Mirny}
\refpages{12--30}
\endreference{}

\reference{shmak17}
\reftype{4}
\refauthf{Shmakov}{I. I.}{}
\refyear{2018}
\reftitle{Problems of scientific support in exploration for diamonds}
\refbook{Geology and minerageny of northern Eurasia. Materials of the
meeting dedicated to the 60-th anniversary of the Institute of Geology
and Geophysics of the Siberian Branch of the Academy of Sciences of the
USSR}
\refpages{265}
\refpubl{Sobolev Institute  of geology and mineralogy SB RAS}
\reflocat{Novosibirsk, Russia}
\endreference{in Russian}

\reference{jur06}
\reftype{3}
\refauthf{Zhuravlev}{Yu. I.}{}
\refauth{Ryazanov}{V. V.}{}
\refauth{Sen'ko}{O. V.}{}
\refyear{2006}
\refbook{Recognition. Mathematical methods. Software system. Practical
applications.}
\refpages{159}
\refpubl{Fazis}
\reflocat{Moscow}
\endreference{in Russian}

\vspace*{1.5\baselineskip}

\end{references}

\end{document}
