AM13 Gaussian, accounting for modeling errors

A Matlab script for the following example is available at sippi_AM13_metropolis_modeling_error.m.

The use of any of the forward models defined above, will be approximation to solving the perfect forward problem. This leads to a 'modeling' error as demonstrated by HCM14. If one has access to an optimal (but perhaps computational inefficient) forward model, and a faster (less accurate) forward model, then a Gaussian model of the modeling error caused by using the approximate, as opposed to the optimal, forward model can be estimated using sippi_compute_modelization_forward_error.

SIPPI allows accounting for such modeling error through the dt and Ct fields for the data structure.

The setup of the data, prior and forward structures is identical to the one described in the previous example.

 %% Load the travel time data set from ARRENAES
 clear all;close all
 D=load('AM13_data.mat');
 options.txt='AM13';

 %% SETUP DATA
 id=1;
 data{id}.d_obs=D.d_obs;
 data{id}.d_std=D.d_std;
 data{id}.Ct=D.Ct+1; % Covariance describing modeling error

 %% SETUP PRIOR
 im=1;
 prior{im}.type='FFTMA';
 prior{im}.name='Velocity (m/ns)';
 prior{im}.m0=0.145;
 prior{im}.Va='.0003 Sph(6)';
 dx=0.15;
 prior{im}.x=[-1:dx:6];
 prior{im}.y=[0:dx:13];
 prior{im}.cax=[.1 .18];

 % SETUP THE FORWARD MODEL USED IN INVERSION
 forward.forward_function='sippi_forward_traveltime';
 forward.sources=D.S;
 forward.receivers=D.R;
 forward.type='fat';forward.linear=1;forward.freq=0.1;

In order to compute the modeling with respect to using the 'Born' type forward model, one can define a new forward structure, here forward_full, and estimate a Gaussian model for the modeling error using

 % SETUP THE 'OPTIMAL' FORWARD MODEL
 forward_full.forward_function='sippi_forward_traveltime';
 forward_full.sources=D.S;
 forward_full.receivers=D.R;
 forward_full.type='Born';forward_full.linear=1;forward_full.freq=0.1;

 % COMPUTE MODELING ERROR DUE TO USE OF forward AS OPPOSED TO forward_full
 N=100;
 [Ct,dt,dd]=sippi_compute_modelization_forward_error(forward_full,forward,prior,data,N);

 % ASSIGN MODELING ERROR TO DATA
 data{1}.dt=dt{1 };
 data{1}.Ct=Ct{1};

Sampling of the posterior can proceed exactly as for the previous example, using

 options.mcmc.nite=500000; % optional, default:nite=30000
 options.mcmc.i_sample=500; % optional, default:i_sample=500;
 options.mcmc.i_plot=1000; % optional, default:i_plot=50;
 options=sippi_metropolis(data,prior,forward,options);

 % plot posterior statistics
 sippi_plot_posterior(options.txt);