The Python library SALib provides several global sensitivity analysis methods (including Morris and Sobol). We will test here if (and how) these methods can be applied to sensitivity analysis of Noddy models.
We will start with the simple two-fault model and then extend the analysis to the Gippsland basin model.
In [1]:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())
Out[1]:
In [1]:
import sys, os
import matplotlib.pyplot as plt
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
os.chdir(r'/Users/flow/git/pynoddy/docs/notebooks/')
repo_path = os.path.realpath('../..')
import pynoddy.history
import pynoddy.events
In [3]:
reload(pynoddy.history)
reload(pynoddy.events)
nm = pynoddy.history.NoddyHistory()
# add stratigraphy
strati_options = {'num_layers' : 8,
'layer_names' : ['layer 1', 'layer 2', 'layer 3', 'layer 4', 'layer 5', 'layer 6', 'layer 7', 'layer 8'],
'layer_thickness' : [1500, 500, 500, 500, 500, 500, 500, 500]}
nm.add_event('stratigraphy', strati_options )
# The following options define the fault geometry:
fault_options = {'name' : 'Fault_W',
'pos' : (4000, 3500, 5000),
'dip_dir' : 90,
'dip' : 60,
'slip' : 1000}
nm.add_event('fault', fault_options)
# The following options define the fault geometry:
fault_options = {'name' : 'Fault_E',
'pos' : (6000, 3500, 5000),
'dip_dir' : 270,
'dip' : 60,
'slip' : 1000}
nm.add_event('fault', fault_options)
history = "two_faults_sensi.his"
nm.write_history_tmp(history)
In [4]:
# compute the model:
output_name = "two_faults_sensi_out"
# Compute the model
pynoddy.compute_model(history, output_name)
In [5]:
# Plot output
reload(pynoddy.output)
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('y', layer_labels = strati_options['layer_names'][::-1],
colorbar = True, title="",
savefig = False)
The NoddyHistory class now has a function to evaluate geological units in a "virtual drillhole", i.e.: a 1-D extraction at a given position (x,y). Basically, a model is created with the same kinematic history, but a very small range (around the position x,y), and a small cube size (default: 1 m) for a "quasi-continuous" model export.
The idea is that the data of this "virtual drillhole" could be used for a meaningful sensitivity analysis, as the data is almost continuous.
Here an example of a 1-D export. First, we re-open the model created above and get some model information:
In [6]:
reload(pynoddy.history)
nm_2 = pynoddy.history.NoddyHistory(history)
nm_2.info()
Let's export the geology at a 1-D profile/ virtual borehole at the centre of the model (x = 5000, y = 3500). Note: the resolution of the exported profile can be defined with an additional keyword (resolution higher than 1 m are possible, as well):
In [7]:
# get 1-D line export at centre of model:
drill_1 = nm_2.get_drillhole_data(5000, 3500, resolution = 1)
A plot of the data looks, so far, pretty boring:
In [8]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(drill_1, np.arange(-len(drill_1),0))
Out[8]:
In [9]:
# set parameter changes in dictionary
# define the elemtents and values of the events to change:
changes_fault_1 = {'Dip' : -20}
changes_fault_2 = {'Dip' : -20}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_2.change_event_params(param_changes)
As the get_drillhole_data method automatically performs the model computation step, we can simply get the updated 1-D profile as:
In [10]:
drill_2 = nm_2.get_drillhole_data(5000, 3500, resolution = 1)
If we now compare the new profile to the previous one we can observe the changes:
In [11]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(drill_1, np.arange(-len(drill_1),0), label='original model')
ax.plot(drill_2, np.arange(-len(drill_1),0), label='changed model')
ax.legend(loc = 'lower left')
Out[11]:
So, as expected, we obtain a different 1-D profile. Important to note here is that the resolution of this 1-D export is a lot higher than the model we would obtain with the full 3-D export so that we can use it as a "quasi-continuous" representation at this one position (x,y). This aspect is important for the sensitivity analysis below.
We can, however, just for the purpose of visualising the changes between these two models, also create a quick comparison of slices through the model:
In [12]:
# compute updated model
updated_history = "changed_fault_model.his"
updated_out = "changed_fault_model_out"
nm_2.write_history(updated_history)
pynoddy.compute_model(updated_history, updated_out)
In [13]:
nout = pynoddy.output.NoddyOutput(output_name)
nout_2 = pynoddy.output.NoddyOutput(updated_out)
fig = plt.figure(figsize=(12,12))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
nout.plot_section('y', layer_labels = strati_options['layer_names'][::-1],
colorbar = False, title="", ax = ax1,
savefig = False)
# plot position of drillhole
ax1.axvline(50, c='gray', lw=2)
nout_2.plot_section('y', layer_labels = strati_options['layer_names'][::-1],
colorbar = False, title="", ax = ax2,
savefig = False)
# plot position of drillhole
ax2.axvline(50, c='gray', lw=2)
Out[13]:
The difference between the two models is also obvious in the cross-section - but the resolution is a lot lower (as can be seen from the stair-step faults). However, exporting the entire model in the same resolution as the 1-D profile is prohibitive for memory and computation time reasons - so we will use the 1-D profiles from now on for uncertainty analysis and sensitivity testing.
The next step is to define an objective function based on the exported 1-D profile. This is a bit tricky to do for a general case. Note, for example, that the profile of the "original model" only contains the layers from 8 to 5, whereas the changed model contains layers 8 to 3!
So, we need a way to describe the "distance" between the two profiles. As a first test, let's use a "stratigraphic distance" where we calculate for each cell the difference between the geological layer id between the changed model and the original model, then sum the absolute values and scale by the number of cells.
We define a distance function:
In [115]:
def distance(profile_ori, profile_changed):
"""Calculate the relative stratigraphic distance along a profile"""
return np.sum(np.abs(profile_ori - profile_changed)) / float(len(profile_ori))
In [15]:
distance(drill_1, drill_2)
Out[15]:
Note: the distance can be larger than 1 (as it should be the case!).
In [16]:
reload(pynoddy.history)
reload(pynoddy)
nm_ori = pynoddy.history.NoddyHistory(history)
And we now use a some more "drillhole" positions:
In [17]:
drill_ori_left = nm_ori.get_drillhole_data(2000, 3500)
drill_ori_fault = nm_ori.get_drillhole_data(4000, 3500)
drill_ori_centre = nm_ori.get_drillhole_data(5000, 3500)
# drill_ori_right = nm_ori.get_drillhole_data(8000, 3500)
In [18]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(drill_ori_left, np.arange(-len(drill_ori_left),0))
ax.plot(drill_ori_fault, np.arange(-len(drill_ori_left),0))
ax.plot(drill_ori_centre, np.arange(-len(drill_ori_left),0))
Out[18]:
We now change the fault dips (to 20 degrees in increments of 1), recopmute the drillhole profiles and calculate the distances:
In [42]:
import copy
nm_changed = copy.deepcopy(nm_ori)
dist_left = [0]
dist_fault = [0]
dist_centre = [0]
for i in range(20):
# define the elemtents and values of the events to change:
changes_fault_1 = {'Dip' : -1}
changes_fault_2 = {'Dip' : -1}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_changed.change_event_params(param_changes)
# Export drillhole profiles
drill_changed_left = nm_changed.get_drillhole_data(2000, 3500)
drill_changed_fault = nm_changed.get_drillhole_data(4000, 3500)
drill_changed_centre = nm_changed.get_drillhole_data(5000, 3500)
# calculate distances and add to arrays
dist_left.append(distance(drill_ori_left, drill_changed_left))
dist_fault.append(distance(drill_ori_fault, drill_changed_fault))
dist_centre.append(distance(drill_ori_centre, drill_changed_centre))
In [44]:
fig = plt.figure(figsize = (8,6))
ax = fig.add_subplot(111)
ax.plot(range(0,21), dist_centre, label = "Centre")
ax.plot(range(0,21),dist_left, label = "Left")
ax.plot(range(0,21),dist_fault, label = "Fault")
ax.legend(loc = "upper left")
ax.set_title("Model distance at three positions for increasing dip")
ax.set_xlabel("Dip change")
ax.set_ylabel("Mean stratigraphic distance")
plt.tight_layout()
Distance is increasing for increasing dip, as expected - and qualitatively, it looks ok, as well: the biggest model differences occur in the range of the fault, then the central Graben, and the least changes on the flanks.
In [46]:
nm_changed = copy.deepcopy(nm_ori)
# define the elemtents and values of the events to change:
changes_fault_1 = {'Dip' : -0.1}
changes_fault_2 = {'Dip' : -0.1}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_changed.change_event_params(param_changes)
# Export drillhole profiles
drill_changed_left = nm_changed.get_drillhole_data(2000, 3500)
drill_changed_fault = nm_changed.get_drillhole_data(4000, 3500)
drill_changed_centre = nm_changed.get_drillhole_data(5000, 3500)
# calculate distances and add to arrays
print distance(drill_ori_left, drill_changed_left)
print distance(drill_ori_fault, drill_changed_fault)
print distance(drill_ori_centre, drill_changed_centre)
A last quick check: let's see how the drillhole resolution affects the distance calculation, i.e.: which resolution do we need to get a stable distance estimate?
We will test three cases: a very small dip change (0.01 degree), a small change (0.1 degree), and a larger change (10 degree). For each of these changes, we will increase the drillhole resolution, starting from 10 m, down to 0.1 m and calculate distances for the drillhole position at the position of the fault:
In [116]:
# use logarithmic steps for resolution:
resolution_inc = np.logspace(0,1,15)
print "Check distance for drillhole resolutions of:"
print resolution_inc
In [117]:
distances = {}
dips = [0.1, 1, 10]
for dip in dips: # initialise dictionary to store results
distances[dip] = []
for res in resolution_inc:
print res
for dip in dips:
# create a copy of the history object
nm_changed = copy.deepcopy(nm_ori)
# define the elemtents and values of the events to change:
changes_fault_1 = {'Dip' : -dip}
changes_fault_2 = {'Dip' : -dip}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_changed.change_event_params(param_changes)
# Export drillhole profiles
# drill_changed_left = nm_changed.get_drillhole_data(2000, 3500)
drill_changed_fault = nm_changed.get_drillhole_data(4000, 3500, resolution = res)
drill_ori_fault = nm_ori.get_drillhole_data(4000, 3500, resolution = res)
# drill_changed_centre = nm_changed.get_drillhole_data(5000, 3500)
# calculate distances and add to arrays
# dist_left.append(distance(drill_ori_left, drill_changed_left))
distances[dip].append(distance(drill_ori_fault, drill_changed_fault))
# dist_centre.append(distance(drill_ori_centre, drill_changed_centre))
In [118]:
print distances
In [119]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.loglog(resolution_inc, distances[10], label = "10$^o$")
ax.loglog(resolution_inc, distances[1], label = "0.1$^o$")
ax.loglog(resolution_inc, distances[0.1], label = "0.01$^o$")
ax.legend()
Out[119]:
Results are a bit odd - to do: check implementation of resolution!
We have got the single parts now for the sensitivity analysis. We are now using the global sensitivity analysis methods of the Python package SALib, available on:
https://github.com/jdherman/SALib
As a start, we will test the sensitivity of the model at each drillhole position separately. As parameters, we will use the parameters of the fault events: dip, dip direction, and slip.
The first step for a sensitivity analysis with SALib is to create a text file with a list of parameters and defined ranges for each parameter. Note that we define changes to the parameter (as this is the way we implement the changes with the NoddyHistory.change_event_params method):
In [126]:
param_file = "params.txt"
params = """dip_change_fault_1 -5.0 5.0
dip_dir_change_fault_1 -5.0 5.0
slip_change_fault_1 -100. 100.
dip_change_fault_2 -5.0 5.0
dip_dir_change_fault_2 -5.0 5.0
slip_change_fault_2 -100. 100."""
f = open(param_file, "w")
f.write(params)
f.close()
Now, we create samples for these parameters. The way the samples are generated depends on the type of sensitivity analysis. We are starting with the Sobol method and create samples:
In [127]:
from SALib.sample import saltelli
In [130]:
param_values = saltelli.sample(10, param_file, calc_second_order = True)
With this command, we create a 2-D array where each row contains a set of parameters for our model.
The number of generated samples depends on (1) the number of input parameters, (2) the type of sensitivity analysis, and (3) the choice of calculating second-order sensitivities (for parameter correlation).
We now use this array and run the parameters through our model. Results are again stored in a text file. For simplicity, we will define a function that takes a parameter array, creates the new model, and calculates the distance to the original model at the position of the drillhole:
In [145]:
# set resolution for analysis:
res = 1 # m
# get original model distance again:
drill_ori_fault = nm_ori.get_drillhole_data(4000, 3500, resolution = res)
def change_model_with_paramter_array(param):
# create a copy of the history object
nm_changed = copy.deepcopy(nm_ori)
# define the elemtents and values of the events to change.
# Note: the order of parameters is as in the text file created above!
changes_fault_1 = {'Dip' : param[0],
'Dip Direction' : param[1],
'Slip' : param[2]}
changes_fault_2 = {'Dip' : param[3],
'Dip Direction' : param[4],
'Slip' : param[5]}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_changed.change_event_params(param_changes)
# Export drillhole profiles
# drill_changed_left = nm_changed.get_drillhole_data(2000, 3500)
drill_changed_fault = nm_changed.get_drillhole_data(4000, 3500, resolution = res)
# drill_changed_centre = nm_changed.get_drillhole_data(5000, 3500)
# calculate and return distance:
return distance(drill_ori_fault, drill_changed_fault)
We can calculate the distance for any parameter set defined in the model:
In [152]:
param_set = 8
dist = change_model_with_paramter_array(param_values[param_set-1])
print("Distance for parameter set %d: %.4f" % (param_set,dist))
Now, we run the entire parameter set through the model and store the results in a text file:
In [154]:
dist_sensi = []
for params in param_values:
dist_sensi.append(change_model_with_paramter_array(params))
Look at sensitivities to check if we actually changed different "quasi-continuous" values:
In [158]:
plot(dist_sensi, 'o-')
Out[158]:
We now save the sensitivities to a file and then perform the sobol sensitivity analysis:
In [162]:
# save results
results_file = 'dist_sensi_graben.txt'
np.savetxt(results_file, dist_sensi, delimiter=' ')
In [163]:
from SALib.analyze import sobol
In [164]:
Si = sobol.analyze(param_file, results_file,
column = 0,
conf_level = 0.95,
calc_second_order = True,
print_to_console=False)
In [217]:
# create composite matrix for sensitivities
n_params = 6
comp_matrix = np.ndarray(shape = (n_params,n_params))
for j in range(n_params):
for i in range(n_params):
if i == j:
comp_matrix[i,j] = Si['S1'][i]
else:
comp_matrix[i,j] = Si['S2'][i,j]
comp_matrix[j,i] = Si['S2'][i,j]
# print comp_matrix
# define labels for figure: phi = dip, d = dip direction, s = slip, subscript = fault event
label_names = ["","$\phi_1$", "$d_1$", "$s_1$", "$\phi_2$", "$d_2$", "$s_2$"]
# Create figure
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.imshow(comp_matrix, interpolation='nearest', cmap='RdBu_r',
vmax = np.max(np.abs(comp_matrix)),
vmin = -np.max(np.abs(comp_matrix)),
)
ax.yaxis.set_ticks_position("both")
ax.xaxis.set_ticks_position("top")
ax.set_xticklabels(label_names)
ax.set_yticklabels(label_names)
# ax.set_title("Sensitivities")
ax.set_xlabel("Parameter Sensitivities")
fig.colorbar(im)
plt.tight_layout()
In [206]:
# re-compute original model, just in case:
ori_history = "fault_ori.his"
ori_out = "fault_ori_out"
nm_ori.write_history(ori_history)
pynoddy.compute_model(ori_history, ori_out)
no_ori = pynoddy.output.NoddyOutput(ori_out)
In [210]:
fig = plt.figure(figsize=(8,6))
ax1 = fig.add_subplot(111)
no_ori.plot_section('y', layer_labels = strati_options['layer_names'][::-1],
colorbar = False, title="", ax = ax1,
savefig = False)
# plot position of drillhole
ax1.axvline(40, c='gray', lw=2);
The position of the "drillhole" is actually at 4000 m (cell 40), coinciding with the position of the first fault:
In [212]:
print("Position of first fault:\t%.2f" % nm_ori.events[2].properties['X'])
print("Position of second fault:\t%.2f" % nm_ori.events[3].properties['X'])
So we could expect that the drillhole geology is sensitive to changes in the parameters of this fault.
Still, the simulated model at this position is more sensitive to the properties of the second fault, at position 6000 (see sensitivity matrix)! This is a bit counter-intuitive, let's see what actually happens: we generate a couple of the realisations for the drawn parameter sets and visualise the cross-sections:
In [214]:
# create base figure:
fig = plt.figure(figsize = (12,16))
for i in range(10):
# create axis for figure
ax = fig.add_subplot(5,2,i+1)
param = param_values[i]
nm_changed = copy.deepcopy(nm_ori)
# read and assign parameter changes:
changes_fault_1 = {'Dip' : param[0],
'Dip Direction' : param[1],
'Slip' : param[2]}
changes_fault_2 = {'Dip' : param[3],
'Dip Direction' : param[4],
'Slip' : param[5]}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_changed.change_event_params(param_changes)
# compute model and create section:
tmp_his = "fault_tmp.his"
tmp_out = "fault_tmp_out"
nm_changed.write_history(tmp_his)
pynoddy.compute_model(tmp_his, tmp_out)
no_changed = pynoddy.output.NoddyOutput(tmp_out)
no_changed.plot_section('y', layer_labels = strati_options['layer_names'][::-1],
colorbar = False, title="", ax = ax,
savefig = False)
# plot position of drillhole
ax.axvline(40, c='gray', lw=2);
If we look at these results, it is quite good to see that the influence of the right (second) fault is more important as it affects the entire lower part of the "drillhole". It is also obvious that the fault dip and the fault slip play a major role (dip direction is actually quite insiginficant as the profile is taken at the center of the model where the fault position is defined).
With these considerations, here the interpretation for the results of the sensitivity analysis:
We analysed the sensitivity of the modelled geology along a 1-D profile (a virtual "drillhole") with respect to fault parameters at two defined fault events. A global sobol sensitivity analysis revealed that:
The result of the sensitivity analysis provides a good insight into the sensitivities of modelled geology along a 1-D profile with respect to event parameters. Results are initially a bit counter-intuitive, but make perfect sense (and so, are actually quite interesting!).
The next step would now be to:
I am actually quite happy with the first test and think that we could get a nice analysis out for the Gippsland Basin - maybe also some interesting counter-intuitive results!
In [ ]: