In this python notebook we explore a facies classification model using a Deep Neural Network arquitecture that takes into account spatial dependencies to outperform the prediction model proposed in the prediction facies from wel logs challenge.
The dataset we will use comes from a class excercise from The University of Kansas on Neural Networks and Fuzzy Systems. This exercise is based on a consortium project to use machine learning techniques to create a reservoir model of the largest gas fields in North America, the Hugoton and Panoma Fields. For more info on the origin of the data, see Bohling and Dubois (2003) and Dubois et al. (2007).
The dataset we will use is log data from nine wells that have been labeled with a facies type based on oberservation of core. We will use this log data to train a classifier to predict facies types.
This data is from the Council Grove gas reservoir in Southwest Kansas. The Panoma Council Grove Field is predominantly a carbonate gas reservoir encompassing 2700 square miles in Southwestern Kansas. This dataset is from nine wells (with 4149 examples), consisting of a set of seven predictor variables and a rock facies (class) for each example vector and validation (test) data (830 examples from two wells) having the same seven predictor variables in the feature vector. Facies are based on examination of cores from nine wells taken vertically at half-foot intervals. Predictor variables include five from wireline log measurements and two geologic constraining variables that are derived from geologic knowledge. These are essentially continuous variables sampled at a half-foot sample rate.
The seven predictor variables are:
The nine discrete facies (classes of rocks) are:
These facies aren't discrete, and gradually blend into one another. Some have neighboring facies that are rather close. Mislabeling within these neighboring facies can be expected to occur. The following table lists the facies, their abbreviated labels and their approximate neighbors.
| Facies | Label | Adjacent Facies |
|---|---|---|
| 1 | SS | 2 |
| 2 | CSiS | 1,3 |
| 3 | FSiS | 2 |
| 4 | SiSh | 5 |
| 5 | MS | 4,6 |
| 6 | WS | 5,7 |
| 7 | D | 6,8 |
| 8 | PS | 6,7,9 |
| 9 | BS | 7,8 |
In [5]:
%%sh
pip install pandas
pip install scikit-learn
pip install keras
In [6]:
from __future__ import print_function
import numpy as np
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
from keras.preprocessing import sequence
from keras.models import Sequential
from keras.layers import Dense, Dropout, Activation
from keras.wrappers.scikit_learn import KerasClassifier
from keras.utils import np_utils
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import KFold
from classification_utilities import display_cm, display_adj_cm
from sklearn.metrics import confusion_matrix, f1_score
from sklearn import preprocessing
We load the training and testing data to preprocess it for further analysis.
In [7]:
filename = 'train_test_data.csv'
training_data = pd.read_csv(filename)
training_data.head(10)
Out[7]:
We fill the missing data values in the PE field with zero and proceed to normalize the data that will be fed into our model.
In [8]:
# Set 'Well Name' and 'Formation' fields as categories
training_data['Well Name'] = training_data['Well Name'].astype('category')
training_data['Formation'] = training_data['Formation'].astype('category')
# Fill missing values and normalize for 'PE' field
training_data['PE'] = training_data['PE'].fillna(value=0)
mean_pe = training_data['PE'].mean()
std_pe = training_data['PE'].std()
training_data['PE'] = (training_data['PE']-mean_pe)/std_pe
# Normalize the rest of fields (GR, ILD_log10, DelthaPHI, PHIND,NM_M,RELPOS)
correct_facies_labels = training_data['Facies'].values
feature_vectors = training_data.drop(['Formation', 'Depth'], axis=1)
well_labels = training_data[['Well Name', 'Facies']].values
data_vectors = feature_vectors.drop(['Well Name', 'Facies'], axis=1).values
scaler = preprocessing.StandardScaler().fit(data_vectors)
scaled_features = scaler.transform(data_vectors)
data_out = np.hstack([well_labels, scaled_features])
In order to start training stage, it is required to format the data by considering a sliding window over the depth component in order to classify a given set of features at some specific depth for each well in the training set.
In [9]:
def preprocess(data_out):
data = data_out
well_data = {}
well_names = list(set(data[:, 0]))
for name in well_names:
well_data[name] = [[], []]
for row in data:
well_data[row[0]][1].append(row[1])
well_data[row[0]][0].append(list(row[2::]))
# Sliding window
positive_lag = 5
negative_lag = 5
chunks = []
chunks_test = []
chunk_length = positive_lag+negative_lag+1
chunks_facies = []
for name in well_names:
if name not in ['STUART', 'CRAWFORD']:
test_well_data = well_data[name]
log_values = np.array(test_well_data[0])
log_values_padded = np.lib.pad(log_values, (negative_lag,positive_lag), 'edge')[:, negative_lag:-positive_lag]
facies_values = np.array(test_well_data[1])
for i in range(log_values.shape[0]):
chunks.append(log_values_padded[i:i+chunk_length, :])
chunks_facies.append(facies_values[i])
else:
test_well_data = well_data[name]
log_values = np.array(test_well_data[0])
log_values_padded = np.lib.pad(log_values, (negative_lag,positive_lag), 'edge')[:, negative_lag:-positive_lag]
facies_values = np.array(test_well_data[1])
for i in range(log_values.shape[0]):
chunks_test.append(log_values_padded[i:i+chunk_length, :])
chunks_facies = np.array(chunks_facies, dtype=np.int32)-1
X_ = np.array(chunks)
X = np.zeros((len(X_),len(X_[0][0]) * len(X_[0])))
for i in range(len(X_)):
X[i,:] = X_[i].flatten()
X_test = np.array(chunks_test)
X_test_out = np.zeros((len(X_test),len(X_test[0][0]) * len(X_test[0])))
for i in range(len(X_test)):
X_test_out[i,:] = X_test[i].flatten()
y = np_utils.to_categorical(chunks_facies)
return X, y, X_test_out
In [10]:
np.random.seed(1337) # For reproducibility
X_train, y_train, X_test = preprocess(data_out)
In [11]:
# Set parameters
input_dim = 77
hidden_dim_1 = 128
hidden_dim_2 = 32
output_dim = 9
batch_size = 32
nb_epoch = 10
In [12]:
def proposed_model():
# Define the model
model = Sequential()
model.add(Dense(128, input_dim=77, init='normal', activation='relu'))
model.add(Dense(32, input_dim=128, init='normal', activation='relu'))
model.add(Dropout(0.3))
model.add(Dense(9, init='normal', activation='softmax'))
# Compile model
model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
return model
In order to evaluate our classification model accurary we will use the our following defined metrics, based on the confusion matrix once the classification is performed. The first metric only considers misclassification error and the second one takes into account the fact that facies could be misclassified if they belong to a same group with similar geological characteristics.
In [13]:
def accuracy(conf):
total_correct = 0.
nb_classes = conf.shape[0]
for i in np.arange(0,nb_classes):
total_correct += conf[i][i]
acc = total_correct/sum(sum(conf))
return acc
adjacent_facies = np.array([[1], [0,2], [1], [4], [3,5], [4,6,7], [5,7], [5,6,8], [6,7]])
facies_labels = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS','WS', 'D','PS', 'BS']
def accuracy_adjacent(conf, adjacent_facies):
nb_classes = conf.shape[0]
total_correct = 0.
for i in np.arange(0,nb_classes):
total_correct += conf[i][i]
for j in adjacent_facies[i]:
total_correct += conf[i][j]
return total_correct / sum(sum(conf))
Once the set of parameters are fixed, the training stage of our model begins. We perform a Cross Validation routine to evaluate the performance of the model.
In [14]:
estimator = KerasClassifier(build_fn=proposed_model, nb_epoch=10, batch_size=50, verbose=0)
# Cross Validation
kfold = KFold(n_splits=5, shuffle=True)
results = cross_val_score(estimator, X_train, y_train, cv=kfold)
print(' Cross Validation Results')
print( results )
In [15]:
plt.hist(results, bins= 5)
plt.title("Model Accuracy")
plt.show()
print ( "Cross Validation Average Accuracy %.3f" % results.mean() )
print ( "Cross Validation Standard Deviation %.3f" % results.std() )
In [16]:
# Load the model
model = proposed_model()
#Train model
model.fit(X_train, y_train, nb_epoch=10, verbose=2, shuffle = True)
Out[16]:
In [17]:
# Predict Values on Training set
y_predicted = model.predict( X_train , batch_size=32, verbose=0)
# Print Report
# Format output [0 - 8 ]
y_ = np.zeros((len(y_train),1))
for i in range(len(y_train)):
y_[i] = np.argmax(y_train[i])
y_predicted_ = np.zeros((len(y_predicted), 1))
for i in range(len(y_predicted)):
y_predicted_[i] = np.argmax( y_predicted[i])
# Confusion Matrix
conf = confusion_matrix(y_, y_predicted_)
# Print Results
print ("\nModel Report")
print ("-Accuracy: %.6f" % ( accuracy(conf) ))
print ("-Adjacent Accuracy: %.6f" % ( accuracy_adjacent(conf, adjacent_facies) ))
print ("\nConfusion Matrix")
display_cm(conf, facies_labels, display_metrics=True, hide_zeros=True)
In [63]:
y_test = model.predict( X_test , batch_size=32, verbose=0)
final_predictions = np.zeros((len(y_test),1))
for i in range(len(y_test)):
final_predictions[i] = np.argmax(y_test[i]) + 1
In [73]:
last_submission = pd.read_csv('Prediction.csv')
ls = last_submission['Facies'].as_matrix()
In [74]:
conf = confusion_matrix(ls, final_predictions)
# Print Results
print ("\nModel Report")
print ("-Accuracy: %.6f" % ( accuracy(conf) ))
print ("-Adjacent Accuracy: %.6f" % ( accuracy_adjacent(conf, adjacent_facies) ))
print ("\nConfusion Matrix")
display_cm(conf, facies_labels, display_metrics=True, hide_zeros=True)
In [77]:
# Store results
test_data = pd.read_csv('../validation_data_nofacies.csv')
test_data['Facies'] = final_predictions
test_data.to_csv('Prediction2.csv')
In [ ]: