Sensitivity Analysis

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:

Evaluate which the most important parameters are, and to
Determine which parameters could, in principle, be inverted with suitable information.

Theory: local sensitivity analysis

Basic considerations:

parameter vector $\vec{p}$
residual vector $\vec{r}$
calculated values at observation points $\vec{z}$
Jacobian matrix $J_{ij} = \frac{\partial \vec{z}}{\partial \vec{p}}$

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.

Defining the responses $\vec{z}$

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)$$

Setting up the base model

For a first test: use simple two-fault model from paper



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)

Define parameter uncertainties

We will start with a sensitivity analysis for the parameters of the fault events.



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









    



 STRATIGRAPHY
 FAULT
 FAULT



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()

Calculate total stratigraphic distance



In [61]:

    
def determine_strati_diff(NO1, NO2):
    """calculate total stratigraphic distance between two models"""
    return np.sum(NO1.block - NO2.block) / float(len(NO1.block))

diff = determine_strati_diff(NO1, NO2)



In [62]:

    
diff









    Out[62]:





-2004.0999999999999

Function to modify parameters



In [56]:

    
# 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 [57]:

    
reload(pynoddy.history)
H2 = pynoddy.history.NoddyHistory(history)
H2.change_event_params(param_changes)









    



 STRATIGRAPHY
 FAULT
 FAULT
{2: {'Dip': -20}, 3: {'Dip': -20}}



In [60]:

    
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()

Full sensitivity analysis

Perform now a full sensitivity analysis for all defined parameters and analyse the output matrix:



In [68]:

    
import copy
new_history = "sensi_tmp.his"
new_output = "sensi_out"
def noddy_sensitivity(his, 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
    changes_list = dict(dict())
    # Step 1:
    # create new parameter list to change model
    for event_id, event_dict in param_changes_file.items(): # iterate over events
        for key, val in his.events[event_id]: # iterate over all properties separately
            param_list.append("event_%d_property_%s" % (event_id, key))
            for i in range(2):
                # calculate positive and negative values
                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))
                if i == 1:
                    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_neg = determine_strati_diff(NO1, NO_tmp)
            # calculate central difference
            central_diff = (dist_pos - dist_neg) / (2. * val)
            distances.append(central_di:setn









    



  File "<ipython-input-68-51a7aa81d055>", line 25
    dist_pos = determine_strati_diff(NO1, NO_tmp))
                                                 ^
SyntaxError: invalid syntax



In [65]:

    
changes_fault_1 = {'Dip' : 10}
changes_fault_2 = {'Dip' : 10}
param_changes = {2 : changes_fault_1,
                 3 : changes_fault_2}



In [66]:

    
print len(param_changes)



In [ ]: