Test here: (local) sensitivity analysis of kinematic parameters with respect to a defined objective function. Aim: test how sensitivity the resulting model is to uncertainties in kinematic parameters to:
Basic considerations:
Numerical estimation of Jacobian matrix with central difference scheme (see Finsterle):
$$J_{ij} = \frac{\partial z_i}{\partial p_j} \approx \frac{z_i(\vec{p}; p_j + \delta p_j) - z_i(\vec{p};p_j - \delta p_j)}{2 \delta p_j}$$where $\delta p_j$ is a small perturbation of parameter $j$, often as a fraction of the value.
A meaningful sensitivity analysis obviously depends on the definition of a suitable response vector $\vec{z}$. Ideally, these responses are related to actual observations. In our case, we first want to determine how sensitive a kinematic structural geological model is with respect to uncertainties in the kinematic parameters. We therefore need calculatable measures that describe variations of the model.
As a first-order assumption, we will use a notation of a stratigraphic distance for discrete subsections of the model, for example in single voxets for the calculated model. We define distance $d$ of a subset $\omega$ as the (discrete) difference between the (discrete) stratigraphic value of an ideal model, $\hat{s}$, to the value of a model realisation $s_i$:
$$d(\omega) = \hat{s} - s_i$$In the first example, we will consider only one response: the overall sum of stratigraphic distances for a model realisation $r$ of all subsets (= voxets, in the practical sense), scaled by the number of subsets (for a subsequent comparison of model discretisations):
$$D_r = \frac{1}{n} \sum_{i=1}^n d(\omega_i)$$Note: mistake before: not considering distances at single nodes but only the sum - this lead to "zero-difference" for simple translation! Now: consider more realistic objective function, squared distance:
$$r = \sqrt{\sum_i (z_{i calc} - z_{i ref})^2}$$
In [14]:
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 [15]:
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 [16]:
output_name = "two_faults_sensi_out"
# Compute the model
pynoddy.compute_model(history, output_name)
In [17]:
# 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)
In [26]:
H1 = pynoddy.history.NoddyHistory(history)
# get the original dip of the fault
dip_ori = H1.events[3].properties['Dip']
# dip_ori1 = H1.events[2].properties['Dip']
# add 10 degrees to dip
add_dip = -20
dip_new = dip_ori + add_dip
# dip_new1 = dip_ori1 + add_dip
# and assign back to properties dictionary:
H1.events[3].properties['Dip'] = dip_new
In [27]:
reload(pynoddy.output)
new_history = "sensi_test_dip_changed.his"
new_output = "sensi_test_dip_changed_out"
H1.write_history(new_history)
pynoddy.compute_model(new_history, new_output)
# load output from both models
NO1 = pynoddy.output.NoddyOutput(output_name)
NO2 = pynoddy.output.NoddyOutput(new_output)
# create basic figure layout
fig = plt.figure(figsize = (15,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
NO1.plot_section('y', position=0, ax = ax1, colorbar=False, title="Dip = %.0f" % dip_ori)
NO2.plot_section('y', position=0, ax = ax2, colorbar=False, title="Dip = %.0f" % dip_new)
plt.show()
In [154]:
# def determine_strati_diff(NO1, NO2):
# """calculate total stratigraphic distance between two models"""
# return np.sum(NO1.block - NO2.block) / float(len(NO1.block))
def determine_strati_diff(NO1, NO2):
"""calculate total stratigraphic distance between two models"""
return np.sqrt(np.sum((NO1.block - NO2.block)**2)) / float(len(NO1.block))
diff = determine_strati_diff(NO1, NO2)
In [155]:
diff
Out[155]:
In [156]:
# set parameter changes in dictionary
changes_fault_1 = {'Dip' : -20}
changes_fault_2 = {'Dip' : -20}
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
In [157]:
reload(pynoddy.history)
H2 = pynoddy.history.NoddyHistory(history)
H2.change_event_params(param_changes)
In [158]:
new_history = "param_dict_changes.his"
new_output = "param_dict_changes_out"
H2.write_history(new_history)
pynoddy.compute_model(new_history, new_output)
# load output from both models
NO1 = pynoddy.output.NoddyOutput(output_name)
NO2 = pynoddy.output.NoddyOutput(new_output)
# create basic figure layout
fig = plt.figure(figsize = (15,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
NO1.plot_section('y', position=0, ax = ax1, colorbar=False, title="Original Model")
NO2.plot_section('y', position=0, ax = ax2, colorbar=False, title="Changed Model")
plt.show()
In [168]:
import copy
new_history = "sensi_tmp.his"
new_output = "sensi_out"
def noddy_sensitivity(history_filename, param_change_vals):
"""Perform noddy sensitivity analysis for a model"""
param_list = [] # list to store parameters for later analysis
distances = [] # list to store calcualted distances
# Step 1:
# create new parameter list to change model
for event_id, event_dict in param_change_vals.items(): # iterate over events
for key, val in event_dict.items(): # iterate over all properties separately
changes_list = dict()
changes_list[event_id] = dict()
param_list.append("event_%d_property_%s" % (event_id, key))
for i in range(2):
# calculate positive and negative values
his = pynoddy.history.NoddyHistory(history_filename)
if i == 0:
changes_list[event_id][key] = val
# set changes
his.change_event_params(changes_list)
# save and calculate model
his.write_history(new_history)
pynoddy.compute_model(new_history, new_output)
# open output and calculate distance
NO_tmp = pynoddy.output.NoddyOutput(new_output)
dist_pos = determine_strati_diff(NO1, NO_tmp)
NO_tmp.plot_section('y', position=0, colorbar=False, title="Dist: %.2f" % dist_pos,
savefig=True, fig_filename="event_%d_property_%s_val_%d.png" % (event_id, key,val))
if i == 1:
changes_list[event_id][key] = -val
his.change_event_params(changes_list)
# save and calculate model
his.write_history(new_history)
pynoddy.compute_model(new_history, new_output)
# open output and calculate distance
NO_tmp = pynoddy.output.NoddyOutput(new_output)
dist_neg = determine_strati_diff(NO1, NO_tmp)
NO_tmp.plot_section('y', position=0, colorbar=False, title="Dist: %.2f" % dist_neg,
savefig=True, fig_filename="event_%d_property_%s_val_%d.png" % (event_id, key,val))
# calculate central difference
central_diff = (dist_pos + dist_neg) / (2.)
distances.append(central_diff)
return param_list, distances
In [191]:
changes_fault_1 = {'Dip' : 1.5,
'Dip Direction' : 10,
'Slip': 100.0,
'X': 500.0}
changes_fault_2 = {'Dip' : 1.5,
'Dip Direction' : 10,
'Slip': 100.0,
'X': 500.0}
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
In [192]:
param_list_1, distances = noddy_sensitivity(history, param_changes)
In [171]:
print param_list_1, distances
In [190]:
d = np.array([distances])
plt.bar(arange(0.6,len(distances),1.), np.array(distances[:]))
Out[190]:
In [185]:
range(len(distances))
Out[185]:
In [173]:
c_zz = d * d.transpose()
In [174]:
imshow(c_zz, interpolation='nearest', cmap='gray_r')
colorbar()
Out[174]:
In [93]:
H_ori.events[2].properties
Out[93]:
In [145]:
N_sensi_out = pynoddy.output.NoddyOutput(new_output)
N_sensi_out.plot_section('y', position=0, colorbar=False, title="Changed Model")
In [ ]: