System and method for resistivity interpretation

[0001] The present invention relates to the characterization of electrical resistivity surveys using artificial intelligence.

[0002] Geophysical surveys are often conducted to determine subsurface properties for various reasons, such as leak detection, analysis of anomalies and forensic structural failure analysis.

[0003] Electrical resistivity tomography (ERT) is a well-known geophysical technique for imaging subsurface structures using electrical resistivity measurements made from a current impressed onto electrodes fixed in the ground.

[0004] Resistivity can be computed if the intensity of a current injected into the ground, and the resulting potential difference established between measurement electrodes are known. Resistivity depends on the nature of the soil, the injected current, the method used to measure the injected current, the resulting potential difference between the electrodes and on the geometry of the electric field created between the electrodes.

[0005] In practice, an array of electrodes is positioned in the ground in a known geometric pattern. The positions of the injector and measurement electrodes is known as “topology”. Many topologies are known for electrode placement in an ERT survey. Examples are the Wenner, Wenner-Schlumberger, pole-dipole and dipole-dipole topologies.

[0006] In each topology, a known current is introduced through a pair of injection electrodes. An electric field results. The resulting electric field is measured as a voltage at measurement electrodes in the array. Resistivity valves at various depths in the ground can then be determined from the known current and the measured voltages.

[0007] Referring to FIG. 1, an example dipole-dipole topology 100 is shown. In the dipole-dipole topology, the spacing between the current injection probes 105 and the voltage measurement probes 107 is the same. The spacing is typically a whole number multiplier “n” of the distance “a”.

[0008] Referring to FIG. 2, an example of a channel of test data in an ERT survey using a dipole-dipole topology 200 will be described. Electrodes 202, 204, 206, 208, 210, 212, 214, 216, 218, 220 and 222 are inserted in ground 224 at an equidistant spacing “a” in a straight line. A known current is injected at current injection probes 230. The resulting differential voltage at electrode pairs 237, 239, 241, 243, 245, 247, 249 and 251, is measured. The apparent resistivity of the subsurface can be determined from these measurements at different depths, and forms a two-dimensional matrix of resistivity values. In this example, impressing a current through injection probes 230 and taking a voltage reading at electrode pairs 237 results in a resistivity value for position 280 at depth 260 and horizontal distance 268. Similarly, a voltage reading at electrode pairs 239 results in a resistivity value for position 281 at depth 261 and horizontal distance 269. Likewise, a voltage reading at electrode pairs 241 results in a resistivity value for position 282 at depth 262 and horizontal distance 270. A voltage reading at electrode pairs 243 results in a resistivity value for position 283 at depth 263 and horizontal distance 271. A voltage reading at electrode pairs 245 results in a resistivity value for position 284 at depth 264 and horizontal distance 272. A voltage reading at electrode pairs 247 results in a resistivity value for position 285 at depth 265 and distance 273. A voltage reading at electrode pairs 249 results in a resistivity value for position 286 at depth 266 and distance 274. A voltage reading at electrode pairs 251 results in a resistivity value for position 287 at depth 267 and distance 275. As can be seen, the farther away the voltage electrodes are from the current injection electrodes, the deeper into the earth the resistivity values are taken. The string of resistivity values for this single placement of the current injection probes can be termed a “channel” of resistivity data.

[0009] Referring to FIG. 3, in order to complete an ERT survey, several (or many) channels of data are required. After each set of voltage measurements is taken, the current injection probes are typically advanced in one direction along the array. For example, from position 302 to positions 304, 306, 308, 310, 312 and 314, in sequence. As different voltage differential measurements are taken, separate channels of resistivity readings are derived and advanced in the same sequence from channel 320, channel 322, to channel 324, to channel 328, to channel 330 and to channel 332. The channels form linear diagonal lines at approximately 45° from an interstitial position between the current injection nodes. The result is a set of differential resistivity values associated with a set of subsurface points directly below the electrodes. The set of differential resistivity values is then subjected to curve smoothing routines which result in a set of continuous equipotential lines on a resistivity “map” that can be interpreted to locate subsurface anomalies.

[0010] Many systems currently exist for carrying out an ERT survey. One such system is described in U.S. Pat. No. 8,019,547 to Bryant. The invention described in the '547 Patent employs reconfigurable nodes located at each electrode in the array operated by a master controller. One novel feature of the invention of the '547 Patent is that the nodes may be set as either current impression nodes or voltage measurement nodes. U.S. Pat. No. 8,019,547 is incorporated herein by reference.

[0011] U.S. Pat. Nos. 8,321,160; 7,788,049; 7,813,883 and 7,386,402 to Bryant similarly disclose the use of a reconfigurable array of electrodes inserted into the ground. U.S. Pat. Nos. 8,321,160; 7,788,049; 7,813,883 and 7,386,402 are incorporated herein by reference.

[0012] U.S. patent application Ser. Nos. 17/302,050; 16/948,059; 16/948,061 and 16/948,062 to Bryant disclose the use of nodes comprised of a hermetically sealed controller housing with controller boards interchangeably connected to electrodes inserted into the ground. U.S. application Ser. Nos. 17/302,050; 16/948,059; 16/948,061 and 16/948,062 are incorporated herein by reference.

[0013] U.S. patent application Ser. No. 17/249,974 to Bryant discloses the use of nodes comprised of a hermetically sealed controller housing with a uniformly radiating electrode geometry, such as a hemispherical geometry. U.S. application Ser. No. 17/249,974 is incorporated herein by reference.

[0014] When an ERT survey is conducted using such systems, a matrix of resistivity values is created and generally stored as a table. A resistivity “model” is then derived using mathematical techniques to minimize error in the data. The resistivity model is then subjected to geostatistical analysis, using kriging or other methods to derive a resistivity map. In general, the resistivity map is a graphical display of the resistivity values modified so that it may be interpreted by a human user to locate and identify subsurface anomalies.

[0015] Referring then to FIGS. 4A and 4B, ERT test survey maps 400 and 401 are representative of a typical resistivity map produced from an ERT survey. In these cases, the topology of each ERT system includes twenty nodes in a dipole-dipole topology and the test were taken at parallel locations. The nodes used to produce map 400 were located at positions 402, 404, 406, 408, 410, 412, 414, 416, 418, 420, 422, 424, 426, 428, 430, and 432. The nodes used to produce map 401 were located at positions 434, 436, 438, 440, 442, 444, 446, 448, 450, 452, 454, 456, 458, 460, 462, and 464. Each of the nodes is placed at internode spacing about 2 meters.

[0016] Notably the maps include equipotential shades indicating varying bands of resistivity at varying depths in the soil. The bands are then interpreted to identify subsurface properties and anomalies, such as moisture, voids, density changes etc. High resistivity typically is interpreted as stone or a void. Low resistivity is typically interpreted as water, or wet soil.

[0017] Currently, resistivity data is manually interpreted based on a resistivity map, such as maps 400 and 401, and the known testing parameters, such array shape, current injected, surface soil properties, typical moisture content, and location. However, the mathematical techniques required to generate the resistivity map are computationally expensive. The high computational expense slows processing time and requires a reasonably high level of bandwidth to compute. Further, even though the resistivity maps may be accurate, it requires many years of geophysical training and experience by a human user to interpret them and even then, they may be misinterpreted. Moreover, interpretation of resistivity data often takes a tremendous amount of human time, which increases exponentially with the size of the ERT survey.

[0018] Thus, there is a need in the art for a means of increasing the accuracy and speed with which resistivity data may be interpreted.

[0040] In the description that follows, like parts are marked throughout the specification and figures with the same numerals, respectively. The figures are not necessarily drawn to scale and may be shown in exaggerated or generalized form in the interest of clarity and conciseness. Unless otherwise specified, all tolerances are plus or minus 20%.

[0041] Referring then to FIG. 5, a preferred embodiment of resistivity interpretation system 500 is shown. System 500 is comprised of administrator device 512, system server 506, and client device 510. System server 506 is operatively connected to database 508. Administrator device 512 and client device 510 are each a smart device, such as a personal computer, tablet, or smart phone.

[0042] Administrator device 512, system server 506, and client device 510 are connected through network 504. Network 504 is a wide area network, such as the internet. Administrator device 512, system server 506, and client device 510 are connected to master controller 514 through network 504. Master controller 514 is further connected to a plurality of nodes, such as node 516, node 518, and node 520. Together, the master controller and the nodes form an ERT survey system, such as those set out in U.S. Pat. Nos. 8,321,160, 7,788,049, 7,813,883, 7,386,402, and 8,019,547 to Bryant. The number of nodes used may vary depending on the size of the geographic region being surveyed and the node topology, as known in the art.

[0043] Resident on system server 506 is graphics generator 522, artificial neural network 524, and convolutional neural network 526. In general, graphics generator 522 consumes unprocessed resistivity data from an ERT survey and creates a resistivity map. Convolutional neural network 526 uses image recognition to generate an interpretation of the subsurface characteristics shown on a resistivity map. Artificial neural network 524 generates an interpretation of subsurface characteristics based on unprocessed resistivity data.

[0044] Referring then to FIGS. 6A and 6B, preferred method 600 for interpreting resistivity data will be described.

[0045] At step 602, an ERT survey test initiation message is generated at master controller 514. At step 604, the test initiation message is transmitted to node 516. At step 606, the test initiation message is transmitted to node 518. At step 608, the test initiation message is transmitted to node 520. A test initiation message is transmitted to all nodes in the array.

[0046] At steps 610, 611, and 612 an ERT survey is conducted at each of the nodes in the array. At steps 614, 616, and 618 the resistivity readings taken during the ERT survey are transmitted to master controller 514. At step 620, the master controller logs the test results.

[0047] At step 622, the ERT test results are forwarded to system server 506. At step 624, the system server stores the ERT test results in the database. At step 626, the system server generates a test parameters matrix. The test parameters matrix contains the known environmental factors and the test settings utilized. Environmental factors include characteristics of the survey location, such as the longitude and latitude of the survey location, soil types, weather history, surface conditions, humidity, the soil moisture content, typical regional climate, month of the year, time of day, current weather conditions, and population. Other environmental factors may be included in the matrix. Test settings include information such as the array topology, the minimum and maximum current and voltage levels, error rate, and data acquisition frequency. Other test settings may be included in the matrix.

[0048] In a preferred embodiment, test parameters are included to provide context to resistivity data. The additional context aides in manual evaluation of the data and likewise reduces the computational expense of evaluating the data. However, it should be appreciated that the test parameters are not required to evaluate the data.

[0049] At step 628, the system generates a resistivity matrix representing the resistivity readings for the test. At step 629, the graphics generator creates a resistivity map based on the resistivity readings. At step 630, the resistivity map is transmitted to administrator device 512. At step 631, the resistivity map is manually analyzed. At step 632, an interpretation matrix is generated. In a preferred embodiment, an interpretation matrix includes a list of potential subsurface characteristics that are identified during manual analysis of the resistivity maps. Subsurface characteristics include high density probability, low density probability, moisture probability, and geological anomaly probability. The subsurface characteristics may include more detailed probabilities such as the probability that the following characteristics are present: tunnels, water leaks, stone, clay, minerals, pipes, aquifers, etc., other subsurface characteristics may be included.

[0050] At step 633, the interpretation matrix is transmitted to system server 506. At step 634, the interpretation matrix is associated with the corresponding ERT resistivity test file and the ERT resistivity map, and stored in the database.

[0051] At step 635, the server runs convolutional neural network 526 to produce an output interpretation matrix, as will be further described. At step 638, the output interpretation matrix is returned to the administrator device.

[0052] At step 640, optionally, the administrator device receives a manual verification whether or not the prediction made by the output interpretation matrix is accurate. At step 642, the manual verification is returned to the system server. At step 644, the server retrains the convolutional neural network based on the verification. The system may initiate a retraining session on command or when a predetermined number of new verified resistivity maps are generated. At step 646, the system updates the convolutional neural network synapse weights, as will be further described.

[0053] At step 648, the server runs artificial neural network 524 to produce an output interpretation matrix, as will be further described. At step 652, the interpretation matrix is returned to the administrator device.

[0054] At step 654, the administrator device optionally receives a manual verification of the interpretation matrix. At step 656, the manual verification is returned to the server. At step 658, the server retrains the artificial neural network based on the verification. In a preferred embodiment, the interpretation matrix and its associated resistivity data set are used to automatically train the artificial neural network. At step 660, the system updates the artificial neural network weights, as will be further described.

[0055] At step 662, the system server receives a new ERT dataset including an ERT resistivity map and associated test parameters matrix. At step 664, the system server inputs the new ERT dataset into and runs trained convolutional neural network 526 to characterize the resistivity data, as will be further described. At step 665, the system runs trained artificial neural network 524 to characterize the resistivity data, as will be further described.

[0056] At step 666, the system server generates a report containing the output interpretation matrix. At step 668, the report is transmitted to client device 510.

[0057] Referring then to FIG. 7A, a preferred embodiment of a convolutional neural network for interpreting subsurface characteristics in a resistivity map will be further described.

[0058] Neural network 700 includes input layer 702, hidden layer 704, and output layer 708. Other numbers of nodes and layers may be used.

[0059] Input layer 702 has one node for each input value (IV) considered in interpreting a resistivity map. For instance, the nodes may represent extracted values from a resistivity map, and test parameters, as will be further described. The input layer may have a greater or lesser number of nodes based on the number of input values. In this illustrative example, a 10×1 input matrix is shown, but more input values are typically used to accurately represent the resistivity map and test parameters, as will be further described. Training data also includes a known output interpretation matrix for each data set, as will be further described.

[0060] Each node of input layer 702 is connected to each node of hidden layer 704 via α synapses. During training, each α synapse is initially assigned a random weight between 0 and 1. The input layer nodes are multiplied by the α synapses, summed, and processed using an activation function to produce a hidden value, HV, for hidden layer nodes, as will be further described.

[0061] Hidden layer 704 is comprised of a plurality of weighted nodes. Hidden layer 704 preferably includes a bias node. The inclusion of a bias node reduces the variance and increases efficiency. In this example, convolutional neural network is shown with eleven (11) nodes in the hidden layer, however, one of skill in the art will recognize that other arrangements, numbers and layers of neurons are possible that may provide the desired predictive features of the invention.

[0062] Each node of hidden layer 704 is connected to output layer 708 by β synapses. During training, each β synapse is initially assigned a random weight between 0 and 1.

[0063] The nodes of hidden layer 704 are multiplied by β synapses, summed, and processed using an activation function to produce the output probability values for output layer 708, as will be further described.

[0064] Neural network 700 includes a plurality of output nodes in output layer 708. Each output node represents a probability that a specific subsurface characteristic, feature, or anomaly, is present in a resistivity map. In this example, output layer 708 is shown with output node 710, output node 712, and output node 714. The values F1, F2, and F3 for each of nodes 710, 712, and 714, respectively, might range from 0 to 1 with a value of 1 being 100% probability that the characteristic is in the survey site. However, it should be appreciated that any number of output nodes could be used to represent a plurality of possible subsurface characteristics.

[0065] Referring to FIG. 7B, a preferred process for training convolutional neural network 700 will be further described.

[0066] Resistivity matrix 716 represents the resistivity data produced in an ERT survey. Resistivity matrix 716 is input into graphics generator 718 resident on the system server to produce resistivity map 720, as previously described.

[0067] Test parameters matrix 722 is comprised of the known environmental factors and the test settings utilized for the ERT survey. Examples of environmental factors include details of the physical location and geology being tested, such as the longitude and latitude coordinates, the soil type, the moisture content of the soil, the weather conditions, time of day etc. Examples of test settings include details of the ERT system being used to output the survey, such as, array type, and minimum and maximum current and voltage, as previously described. In a preferred embodiment, all test parameters are provided on a normalized scale from 0 to 1. Non-numerical parameters, such as array type and soil type, are assigned numerical values from a preferred table. For example, possible values representing array types are listed in Table 1 below.

[0068]

	TABLE 1
		Test
		Parameters
		Matrix
	Array	Value
	Dipole-Dipole Array	0.1
	Pole-Dipole Array	0.2
	Wenner Array	0.3
	Schlumberger Array	0.4

[0069] As another example, soil type may be represented by the percentage of each soil type in the soil composition found at the ERT survey site, exemplary values are listed in Table 2 below. In this example, the percentages, when totaled, should equal 100%.

[0070]

	TABLE 2
		Test
		Parameters
		Matrix
	Soil Type	Value
	Sand	0.05
	Loam	0.84
	Clay	0.08
	Rock	0.03

[0071] Interpretation matrix 724 is a list of probabilities that each of a set of subsurface characteristics, such as tunnels, water leak, stone, clay, minerals, oil, aquifers, etc., is present in the survey site represented by the resistivity map. In one embodiment, the interpretation matrix is manually generated based on actual knowledge of whether or not a characteristic is present. In this embodiment, the characteristic is either 100% present or 0% present based on actual samples taken and borings conducted before or after the ERT survey. In another embodiment, the interpretation matrix is a list of probabilities ranging between 0% and 100% generated based on statistical analysis and expert human interpretations of resistivity data. A combination of both scenarios may also be used. In either case, the values used for interpretation matrix 724 are considered known values. Exemplary output values and associated characteristics are shown in Table 3 below.

[0072]

	TABLE 3
		Interpretation
		Matrix
	Characteristic	Value
	Tunnel	20%
	Water Leak	5%
	Limestone	100%
	Oil	0%
	. . .	. . .

[0073] During training, resistivity map 720, test parameters matrix 722, and interpretation matrix 724, are concatenated and input into neural network 700 to derive a set of trained α and β synapse weights and training output 726.

[0074] During training, the known probability values in interpretation matrix 724 are compared to the derived probability values in output layer 708 to determine a margin of error. During backpropagation, the errors attributable to the α and β synapse weights are determined and the weights are adjusted, as will be further described. The convolutional neural network generates a new interpretation matrix using the new set of weights. This process is repeated until an optimal loss is achieved to minimize the margin of error. When optimal loss is achieved, the set of ideal trained α and β synapse weights, training output 726, is stored in the system server, as will be further described. Training output 726 is used in the functioning convolutional neural network, as will be further described.

[0075] Referring then to FIG. 7C, a process of using a trained convolutional neural network 701 will be described.

[0076] Resistivity matrix 717 is input into graphics generator 719 resident on the system server to produce resistivity map 721, as previously described.

[0077] Resistivity map 721 and test parameters matrix 723 are concatenated and input into trained convolutional neural network 701 to derive output interpretation matrix 728.

[0078] Referring then to FIGS. 7D and 7E, feature extraction method 760 for processing images for input into convolutional neural network 700 will be further described.

[0079] The processes used in feature extraction method 760 include image input layer 730, convolution layer 732, and pooling layer 734. The processes systematically extract information from the resistivity maps and format the data for input layer 702 of convolutional neural network 700.

[0080] Image input layer 730 is preferably comprised of a large data set of ordered pixels which make up a resistivity map image. Each resistivity map image is associated with one test parameters matrix. In a preferred embodiment, the input resistivity map images are the same size and format. In this illustrative example, input image 738 is a 5×5 pixel image. One of skill in the art will recognize that typical resistivity map images are comprised of a much larger pixel dimensions in order to characterize resistivity data, but convolutional neural network 700 still processes images in the same manner as explained in this example.

[0081] At step 762, each convolutional layer 732 scans the resistivity map for specific graphical features and performs a dot matrix multiplication to simplify the data to create a feature map. There may be many convolutional layers and many feature maps. Each convolution layer scans for different features in the image, such as edges, patterns, gradients, red, green, and blue, to produce a corresponding feature map, such as feature map 746. Features may include specific anomalies recognized in resistivity maps, such as tunnels and leaks. It should be appreciated that more convolution layers and feature maps may be utilized with increasing specificity in feature extraction.

[0082] In this example, feature map 746 is a 3×3 pixel matrix created by scanning regions of input image 738 using feature window 736. A feature window is defined as a set number of pixels per region used to examine an image. In this example, feature window 736 is configured to determine whether or not the section of the image examined has a single value, such as the color red, represented as the digit 1. Feature window 736 is configured as follows:

[0083]

1	1	1
1	1	1
1	1	1

[0084] Different feature windows may be configured in alternate ways to look for various features. For instance, a feature window searching for a vertical line may be configured as follows:

[0085]

0	1	0
0	1	0
0	1	0

[0086] Feature window 736 is overlayed on region 737 of image 738 in a bitwise multiplication. In region 737, red is represented by 1 and non-red is 0. Region 737 is represented below:

[0087]

1	0	1
0	1	0
0	1	1

[0088] Feature window 736 moves across image 738 by defined horizontal stride 740 and vertical stride 741, which are both one pixel. Other stride dimensions may be used. The stride is the unit length that the window is moved horizontally and vertically for each region. The convolution layer scans the entire image by region, moving window 736 horizontally by stride 740 until all pixel columns are considered and then advancing downward by stride 741 to repeat scanning horizontally until all columns and rows are considered. As the window moves across the image, a matrix multiplication of the window and the image is computed and stored in the appropriate region position of the feature map.

[0089] For instance, feature window 736 examines region 737 for the color red to produce region 742 of feature map 746. Region 742 is calculated as follows:
(1×1)+(0×1)+(1×1)+(0×1)+(1×1)+(0×1)+(0×1)+(1×1)+(1×1)=5

[0090] At step 764, pooling layer 734 is used to compress the spatial information for each feature map and reduce overfitting the network. The feature maps are compressed in the same order that they were created. In this example, pooling layer 734 scans feature map 746 using exemplary window 744 to produce pooled feature map 750. Window 744 is a 2×2 matrix window. The window moves across feature map 746 horizontally by a stride of one and then vertically by a stride of one until all regions are considered. In this example, a “stride” includes a single column width and row height. In this example, window 744 is configured to compress the values of the feature map into four regions using max pooling. In “max pooling” the maximum value in the window is located and transferred to the corresponding region of the pooled feature map. In this example, region 748 results from window 744 and is “5”. In an alternate embodiment, “average pooling” may be used, which takes the average value in a region.

[0091] In yet another embodiment, “global feature mapping” is utilized to produce uniform compressed feature maps in pooling layer 734. By using global feature mapping, convolutional neural network 700 is able to accept images of different sizes, such as with various array sizes. However, global feature mapping is more computationally expensive. Thus, in a preferred embodiment, images are reduced to a standard size before being input into the neural network to increase efficiency. In a preferred embodiment, the images are reduced to the same size but still contains images representing resistivity surveys conducted using various array sizes and data representing the array size and configuration are input into a fully connected input layer. In another embodiment, a different neural network may be employed for each different size or configuration of resistivity maps in order to standardize the number of input values.

[0092] It will be appreciated by those skilled in the art that other configurations of window size, stride size, number of feature maps, convolutional layers, and pooling layers may be employed to gain the advantages of the invention.

[0093] At step 766, the pooled data is flatted into a single column matrix. “Flattening” the data arranges the rows and columns into a single column matrix. Pooled feature maps are flattened in the same order as their corresponding feature maps are created. Each pooled feature map is parsed out according to a set directionality, generally left to right top to bottom. In this example, the values from pooled feature map 750 are flattened into a 4×1 matrix from left to right top to bottom. So in this example, the values from six (6) maps in pooling layer 734 would be flattened into a 24×1 matrix representing of all the extracted features present in the image.

[0094] In a preferred embodiment, the values for the test parameters matrix associated with a single resistivity map are appended to the end of the flattened data to form a single column matrix to completely populate the input layer. The data is then fed forward through convolutional neural network 700, as will be further described.

[0095] In a preferred embodiment, the system uses transfer learning to accelerate training the neural network. Transfer learning utilizes a pretrained convolutional neural network, such as MobileNet V2 trained on the MNIST database, for feature extraction method 760. The benefit of transfer learning is that a model trained on a large dataset can skip learning standard features needed for image recognition, such as numbers, color, edges, patterns, and gradients, and additional features specific to the purpose are integrated and trained, such as resistivity, tunnels, and leaks. The network can then be trained on a large dataset of resistivity maps linked with associated test parameters.

[0096] It should be appreciated that different configurations of the number of convolution and pooling layers may be used. Similarly, any configuration of hyperparameters such as the number of hidden layers, the number of neurons per hidden layer, and type of activation function may be used to gain the advantages of the invention.

[0097] Referring then to FIGS. 7, 8A, 8B, and 8C, a flow chart of method 800 for training convolutional neural network 700 will be described.

[0098] At step 802, optionally, resistivity maps are standardized to a specified size and format.

[0099] At step 804, the images are stored in a specified dataset for the convolutional neural network. In a preferred embodiment, 10,000 to 20,000 images are required to approach an optimal set of synapse values.

[0100] At step 806, features are extracted from a resistivity image and the data is flattened according to feature extraction method 760, as previously described.

[0101] At step 808, a test parameters matrix and an interpretation matrix are associated with each image file.

[0102] Example test parameters matrices and interpretation matrices are shown in Table 4. For this example, the probability values are input on a scale of 0-100% that a subsurface characteristic is present in Resistivity Map ‘X’. Each test parameter matrix and interpretation matrix are associated with a single resistivity map.

[0103]

TABLE 4
Resis-
tivity
Map
(x)	Test Parameters Matrix (y)	Interpretation Matrix (z)
1	[21, 6, 7.5, 30, 150, . . . ]	[10, 45, 75, 22, 85, 0 . . . ]
2	[76, 500, 350, 36, 220 . . . ]	[13, 23, 90, 62, 33, 24 . . . ]
3	[21, 36, 50, 70, 10, . . . ]	[86, 35, 68, 29, 52, 12 . . . ]
. . .	. . .	. . .
x	[y1n, y2n, y3n, y4n, y5n,	[z1n, z2n, z3n, z4n, z5n, z5n,
	. . . ]	. . . ]

[0104] At step 812, the test parameters matrices and the interpretation matrices are normalized to a value between 0 and 1 by dividing each variable by the maximum permitted value of the variables.

[0105] At step 814, the flattened image data is combined with the normalized test parameters to form a fully connected input layer. In a preferred embodiment, for each flattened image file the associated test parameters column matrix is appended to the end of the flattened data.

[0106] At step 828, each of the α synapses are assigned a random weight, α_n, between 0 and 1. At step 830, each of the β synapses are assigned a random weight, β_n, between 0 and 1.

[0107] At step 832, the input layer values and the α synapse weights are multiplied in a matrix operation and summed to determine a first set of hidden values, μ, for each node n in hidden layer 704, according to the following equation:
μ_mn=ΣIV_n×α_m  Eq. 1

• • [0108] where:
    • [0109] n=number of input layer nodes; and
    • [0110] m=number of hidden nodes.

[0111] At step 834, the activation function is applied to the first set of hidden values, μ, to derive the final hidden value of the node, HV. In a preferred embodiment, the activation function is the Sigmoid function. The Sigmoid function is preferred because its derivative can be efficiently calculated. The activation function is shown below:

[0112] $\begin{array}{cc}H{V}_n=S\left({\mu }_{mn}\right)=\frac{1}{1+e^{-{\mu }_{mn}}}=\frac{e^{{\mu }_{mn}}}{e^{{\mu }_{mn}}+1}& \mathrm{Eq}.2\end{array}$

• • [0113] where:
    • [0114] μ_n is the first value of the nodes in the hidden layer; and
    • [0115] HV_n is the final hidden value of the hidden layer nodes.

[0116] At step 836, the values HV_n and the β synapse weights are multiplied in a matrix operation and summed to determine the first value, γ, for the nodes in output layer 708, according to the following equation:
γ=ΣHV_n×β_m  Eq. 3

[0117] At step 838, the activation function is applied to the value, γ, calculated in step 836 to produce the output values for output layer 708. The output values are the predicted probability that a subsurface characteristic is present, as previously described.

[0118] In steps 840-852, during training the neural network backpropagates to minimize Error, as will be further described.

[0119] At step 840, the error in the output is calculated using a mean sum squared loss function. The mean sum squared loss function is the sum for all data points of the square of the difference between the predicted and actual target values divided by the number of sets of data, according to the following equation:

[0120] $\begin{array}{cc}\mathrm{Error}=\sum \frac{{\left(\mathrm{Outpu}{t}_n-z_n\right)}^2}{n}& \mathrm{Eq}.4\end{array}$

• • [0121] where:
    • [0122] Output_n is the predicted probability of a subsurface characteristic; and
    • [0123] Z_n is the known probability value input in the system in step 814.

[0124] At step 842, the margin of error of output layer 708, error margin, is calculated according to the following equation:
Error Margin=y_n−Output_n  Eq. 5

[0125] At step 844, the delta output sum for output layer 708, Δ₁, is calculated. The delta output sum, Δ₁, is calculated by applying the derivative of the sigmoid activation function to the output error calculated in step 840. The derivative of the sigmoid activation function for an output of “x” is x(x−1).

[0126] At step 846, the error attributable to hidden layer 704, the hidden error, is determined. The hidden error is calculated by calculating the dot product of the delta output sum, Δ₁, and the β synapse weights.

[0127] At step 848, the delta output sum for hidden layer 704, Δ₂, is calculated by applying the derivative of the sigmoid activation function to the hidden error calculated in step 846, as previously described.

[0128] At step 850, the a synapse weights are adjusted by calculating the dot product of input layer 702 nodes and the delta output sum, Δ₂. At step 852, the β synapse weights are adjusted by calculating the dot product of hidden layer 704 nodes and the delta output sum, Δ₁. The adjusted weights are returned, and the neural network uses the new weights in the next training iteration. When the network has trained the preset number of iterations the weights are fixed at the values with the lowest mean sum squared loss error.

[0129] At step 854, an automatic evaluation of the trained model is performed and plotted to determine the overall performance of the trained model.

[0130] Steps 832 through 854 are repeated for a preset number of iterations until the convolutional neural network is trained. In a preferred embodiment, a preset number of iterations is used, anywhere from 10,000 to 200,000. Once the system executes the number of iterations, the neural network is considered “trained” and the ideal values of α and β synapses are stored. In an alternate embodiment, if an ideal error, such as 0.01%, is achieved prior to executing all iterations, the neural network is similarly considered trained. Other iterations counts and ideal errors may be used. A higher iteration count reduces the Error and increases the accuracy of the synapse weights.

[0131] At step 858, the synapse weights are stored in the database when training is complete. At step 860, training ends and the convolutional neural network is updated.

[0132] Referring then to FIG. 8D, method 801 for running trained convolutional neural network 701 will be described.

[0133] At step 862, a new ERT survey is conducted and raw resistivity data is transmitted to the system server, as previously described.

[0134] At step 864, a new test parameters matrix is generated, as previously described. At step 866, the resistivity data is processed through the graphics generator to generate a resistivity map, as previously described.

[0135] At step 872, the resistivity map is processed through feature extraction method 760, as previously described. At step 874, the flattened data extracted in step 872 is combined with the test parameters to form a fully connected input layer. In a preferred embodiment, the test parameters matrix values are appended to the end of the flattened data as a normalized column matrix.

[0136] At step 876, the input layer values and the trained α synapse weights are multiplied in a matrix operation and summed to determine a first set of hidden values, μ, for each node n in hidden layer 704, according to equation 1, as previously described.

[0137] At step 878, the activation function is applied to the first set of hidden values, μ, to derive the final hidden value of the node, HV. In a preferred embodiment, the activation function is the Sigmoid function shown in equation 2.

[0138] At step 880, the values HV_n and the trained β synapse weights are multiplied in a matrix operation and summed to determine the first value, γ, for the nodes in output layer 708, according to equation 3, as previously described.

[0139] At step 882, the activation function is applied to the value, γ, calculated in step 880 to produce the output values for output layer 708. The output values are the predicted probability that a subsurface characteristic is present, as previously described.

[0140] At step 884, the system server generates output interpretation matrix 728 using the output values calculated in step 882. The output interpretation matrix list each characteristic and the associated output probability value, that is, the probability that the specified characteristic is present in the resistivity map.

[0141] At step 886, the output interpretation matrix is stored in the system server and returned to the client device or administrator device.

[0142] An example of computer code written in Python to perform one example of the method is shown in FIGS. 9A and 9B. Of course, other code may be used to implement this and other embodiments of the convolutional neural network described.

[0143] Referring then to FIG. 9A, at lines 1-11, the libraries used to run the convolutional neural network are imported. At line 2, Tensorflow is imported to train the neural network. In this example, at line 3, the MNIST dataset is imported from the Keras library to train the convolutional neural network using transfer learning. Additional or different datasets and libraries may be imported and used to achieve similar results. At line 4, the PIL method is called to import the user created dataset into the network.

[0144] At lines 5-10, the Keras methods used to build the convolutional network are imported. The method Dense includes the methods for activation functions. The Conv2D method builds the convolutional layers which search for specific features in an image. The Flatten method flattens the extracted values from the image into a single column matrix. MaxPooling2D reduces the special dimensions of the extracted data.

[0145] At line 11, the NumPy library is imported to perform advanced mathematical functions, such as matrix multiplication, which is needed for neural networks.

[0146] At lines 14-15, a random seed is setup to control the starting point of the randomly assigned weight values to control the evolution of the neural network.

[0147] At line 19, the number of subsurface characteristics, interpretations, are defined. In this example, the number of subsurface characteristics is listed as 10, though the system can be created to define more or less subsurface characteristics.

[0148] At line 20, the image size is defined as 28×28 pixels. This is only an example, a preferred resistivity map image size is 250×1500 pixels.

[0149] At lines 24-26, the datasets are imported and split into train and test sets. At lines 25-26, the images are optionally resized.

[0150] At lines 31-32, the image set is converted to a float data type. At lines 35-36, the images are rescaled, or normalized between 0 and 1.

[0151] At lines 39-40, the test parameters are returned in a binary matrix.

[0152] At lines 42-43, the interpretation matrix is returned in a binary matrix.

[0153] Referring then to FIG. 9B, at line 48, a first convolutional layer is defined. In this example the convolutional layer is set to filter for 32 features in an image, using a 3×3 pixel window, a stride of 1, and utilizes the sigmoid function to normalize the data. However, different numbers and sizes of filters, and other activation functions may be utilized. At line 50, the max pooling method is defined as a 2×2 matrix.

[0154] At lines 53-55, a secondary convolutional layer and pooling layer are defined. In this example, the second convolutional layer is set with a larger number of filters. More or less convolutional layers and pooling layers may be utilized, as previously described.

[0155] At lines 57-60, the data is flattened and normalized using a sigmoid activation function. At line 63, the method Dropout is used with a rate of 0.5 to discard datapoints in order to reduce overfitting the convolutional neural network.

[0156] At line 66, the system converts the output into a probability that a subsurface characteristic is shown in the image, the interpretation matrix. At line 69, the system is compiled and the performance is measured.

[0157] At lines 72-73, the model is trained for 150,000 epochs. At line 75, the loss is calculated as a score for the performance of the convolutional neural network.

[0158] At lines 78-80, the test performance is printed for review.

[0159] At lines 83-86, the neural network is run and the predicted interpretation matrix and the input test parameters matrix are printed for review.

[0160] At line 83, the neural network is run and results are printed. At line 86, actual test results are printed.

[0161] Referring then to FIG. 10A, a preferred embodiment of a process for training an artificial neural network for interpreting subsurface characteristics from raw resistivity data and test parameters will be further described. Artificial neural network 1000 includes input layer 1020, hidden layer 1022, and output layer 1024. Other numbers of nodes and layers may be used, as previously described.

[0162] Input layer 1020 has one node for each input value (IV) considered in interpreting a resistivity data. The input values considered in input layer 1020 are obtained from resistivity matrix 1004 and test parameters matrix 1008, as will be further described.

[0163] Each node of input layer 1020 is connected to each node of hidden layer 1022 via α synapses having an assigned weight between 0 and 1. The input layer nodes are multiplied by the α synapses, summed, and processed using an activation function to produce a hidden value, HV, for hidden layer nodes, as will be further described.

[0164] Hidden layer 1022 is comprised of a plurality of weighted nodes. Hidden layer 1022 preferably includes a bias node. Each node of hidden layer 1022 is connected to output layer 1024 by β synapses assigned a weight between 0 and 1. The nodes of hidden layer 1022 are multiplied by β synapses, summed, and processed using an activation function to produce the output probability values for output layer 1024, as will be further described.

[0165] In this example, output layer 1024 is shown with three (3) nodes 1012, 1016, and 1018, representing the probabilities that three different subsurface characteristics are present.

[0166] During training, resistivity matrix 1004, test parameters matrix 1008, and interpretation matrix 1010, are input into artificial neural network 1000 to derive of trained α and β synapse weights, training output 1014, as will be further described.

[0167] In a preferred embodiment, resistivity matrix 1004 is input into graphics generator 1002 resident on the system server to produce resistivity map 1006. In one embodiment, resistivity map 1006 is processed through the convolutional neural network to derive interpretation matrix 1010 to train neural network 1000. In another embodiment, the resistivity map is manually interpreted to derive interpretation matrix 1010.

[0168] During training, the known probability values in interpretation matrix 1010 are compared to the derived probability values in output layer 1024 to determine a margin of error and the α and β synapse weights are adjusted through backpropagation, as will be further described. When optimal loss is achieved, the set of ideal trained α and β synapse weights, training output 1014, is stored in the system server and used in the functioning artificial neural network, as will be further described.

[0169] Referring then to FIG. 10B, a process of use of trained artificial neural network 1001 will be described.

[0170] Resistivity matrix 1005 and test parameters matrix 1009 are input into trained artificial neural network 1001 to derive interpretation matrix 1026. In one embodiment, trained neural network 1001 includes input layer 1020. Input layer 1020 is operatively connected to hidden layer 1022 by synapses α. Hidden layer 1022 is connected to output layer 1024 by synapses β. The output layer produces output interpretation matrix 1026. The output interpretation matrix provides a set of probability values that each specified characteristic, such as a leak or a tunnel, is present in the subsurface geology that produced the resistivity matrix. Notably, the necessity to generate a resistivity map is eliminated in this embodiment, drastically reducing computational cost required and increasing the speed at which an interpretation of subsurface features can be made.

[0171] Referring then to FIGS. 11A and 11B, a flow chart of method 1100 for training artificial neural network 1000 will be described.

[0172] At step 1110, training data is input into the artificial neural network. The training data is supplied to the algorithm as sets of single column matrices representing the resistivity matrix, the test parameters matrix, and the interpretation matrix, as previously described. In a preferred embodiment, in each dataset, the test parameters matrix is appended to the end of the resistivity matrix to form a single n×1 matrix.

[0173] At step 1112, in a preferred embodiment, all of the values input in step 1110 are normalized to a value between 0 and 1 by dividing each variable by the maximum permitted value of the variables.

[0174] At step 1114, each of the α synapses are assigned a random weight, α_n, between 0 and 1. At step 1116, each of the β synapses are assigned a random weight, β_n, between 0 and 1.

[0175] At step 1118, the input layer values and the α synapse weights are multiplied in a matrix operation and summed to determine a first set of hidden values, μ, for each node n in hidden layer 1022 according to equation 1, as previously described.

[0176] At step 1120, the sigmoid activation function is applied to the first set of hidden values, μ, to derive the final hidden value of the node, HV, using equation 2, as previously described.

[0177] At step 1122, the values HV_n and the β synapse weights are multiplied in a matrix operation and summed to determine the first value, γ, for the nodes in output layer 1024, according to equation 3, as previously described.

[0178] At step 1124, the activation function is applied to the value, γ, calculated in step 1122 to produce the output values for output layer 1024. The output values are the predicted probability that a subsurface characteristic is present, as previously described.

[0179] In steps 1126-1138, during training the neural network backpropagates to minimize Error, as will be further described.

[0180] At step 1126, the error in the output is calculated using a mean sum squared loss function according to equation 4, as previously described.

[0181] At step 1128, the margin of error of output layer 1024 is calculated according to equation 5, as previously described.

[0182] At step 1130, the delta output sum for output layer 1024, Δ₁, is calculated by applying the derivative of the sigmoid activation function to the output error calculated in step 1126, as previously described.

[0183] At step 1132, the error attributable to hidden layer 1022, the hidden error, is determined. The hidden error is calculated by calculating the dot product of the delta output sum, Δ₁, and the β synapse weights.

[0184] At step 1134, the delta output sum for hidden layer 1022, Δ₂, is calculated by applying the derivative of the sigmoid activation function to the hidden error calculated in step 1132, as previously described.

[0185] At step 1136, the a synapses weights are adjusted by calculating the dot product of input layer 1020 nodes and the delta output sum, Δ₂. At step 1138, the β synapses weights are adjusted by calculating the dot product of hidden layer 1022 nodes and the delta output sum, Δ₁. The adjusted weights are returned, and the neural network uses the new weights in the next training iteration. When the network has trained the preset number of iterations, the weights are fixed at the values with the lowest mean sum squared loss error.

[0186] At step 1140, an automatic evaluation of the trained model is performed and plotted to determine the overall performance of the trained model.

[0187] Steps 1118 through 1140 are repeated for a preset number of iterations or until an acceptable margin of error is reached.

[0188] At step 1144, the ideal values of the synapse weights are stored in the database when training is complete. At step 1146, training ends and the artificial neural network is updated.

[0189] Referring then to FIG. 11C, method 1101 for using trained artificial neural network 1001 will be described.

[0190] At step 1150, a new ERT survey is conducted and transmitted to the system server, as previously described.

[0191] At step 1151, the resistivity data is normalized to a scale of between 0 and 1.

[0192] At step 1152, a new test parameters matrix is generated. In a preferred embodiment, the test parameters are provided on a 0 to 1 scale, as previously described.

[0193] At step 1154, the resistivity matrix is appended to the test parameters matrix to form a dataset.

[0194] At step 1156, the input layer values and the trained α synapse weights are multiplied in a matrix operation and summed to determine a first set of hidden values, μ, for each node n in hidden layer 1022, according equation 1.

[0195] At step 1158, the activation function is applied to the first set of hidden values, μ, to derive the final hidden value of the node, HV. In a preferred embodiment, the activation function is the Sigmoid function shown in equation 2.

[0196] At step 1160, the values HV_n and the trained β synapse weights are multiplied in a matrix operation and summed to determine the first value, γ, for the nodes in output layer 1024, according to the equation 3.

[0197] At step 1162, the activation function is applied to the value, γ, calculated in step 1160 to produce the output values for output layer 1024. The output values are the predicted probabilities that specific subsurface characteristics are present, as previously described.

[0198] At step 1164, the system server generates output interpretation matrix 1026 using the output values calculated in step 1162. The output interpretation matrix list each characteristic and the associated output probability value.

[0199] At step 1166, the output interpretation matrix is stored in the system server and returned to the client device or administrator device.

[0200] An example of computer code written in Python to perform one example of the method is shown in FIGS. 12A and 12B. Of course, other code may be used to implement this and other embodiments of the artificial neural network described.

[0201] Referring then to FIG. 12A, at line 1, the NumPy library is imported to perform advanced mathematical functions, such as matrix multiplication, which is needed for neural networks. Additional or different libraries may be imported and used, such as Tensorflow and Keras.

[0202] At lines 3-7, the training data is input as a matrix. As an example, four (4) sets of resistivity matrix values and test parameter values, x, are input and three (3) sets of interpretation probability values, y are included. In practice between 1,000 and 25,000 data sets would be employed to train each neural network.

[0203] At lines 9-10, the data is scaled by dividing the x array by its maximum value and they array by 100, or the maximum possible probability value input. At lines 13-15, the data is split into training and testing data. In this example, the first three sets of data are used to train the neural network. Once training is complete, the neural network is tested by deriving an interpretation matrix for the fourth set of resistivity matrix and test parameter data.

[0204] At lines 17-23, the neural network is created. In this example, the network has one input layer with ten (10) nodes, one output layer with eight (8) nodes, and one hidden layer with eleven (11) nodes.

[0205] At lines 25-27, the α synapse and β synapse weights are defined as 11×10 and 11×8 arrays, respectively, and assigned random values.

[0206] At lines 29-35, a forward propagation function is defined. When called, the forward propagation function executes steps 1118 through 1124, as previously described. At lines 36-38, the sigmoid activation function is defined. At lines 40-42, the derivative of the sigmoid activation function is defined.

[0207] Referring then to FIG. 12B, at lines 44-53, the backpropagation function is defined. When called, the back propagation function executes steps 1126 through 1138, as previously described.

[0208] At lines 55-58, the system is set to train by using forward propagation to produce an output. At lines 60-63, the system is set to save the synapse weights. At lines 65-69, the system is set to print the interpretation matrix input and the resulting predicted probability of a subsurface characteristic being present after the neural network is trained.

[0209] At lines 71-84, the system runs the neural network until it is trained. In this example the system will run 150,000 iterations of training unless an ideal loss of 0.001% is achieved. In this example, during training, the neural network is set to print the scaled input, actual output, predicted output, and mean sum squared loss for each iteration run. Once the network is trained, at line 86, the system saves the ideal synapse weights.

[0210] At line 87, the system returns an interpretation matrix for the fourth set of input data.