From b02d68b51ebc4c5455548102281e82fe1f57e4ef Mon Sep 17 00:00:00 2001 From: FreekvanLeijen Date: Tue, 25 Jul 2017 16:44:19 +0200 Subject: [PATCH] Added step3 with a modeling example and updated README file. --- example/README | 41 ++-- .../cupido_example_step3_example_modeling.m | 184 ++++++++++++++++++ 2 files changed, 205 insertions(+), 20 deletions(-) create mode 100644 example/cupido_example_step3_example_modeling.m diff --git a/example/README b/example/README index d721c81..e0b1259 100644 --- a/example/README +++ b/example/README @@ -1,20 +1,21 @@ -CUPiDO - Connecting Undifferenced Points in Deformation Observations - -This directory contains an example setup to demonstrate the CUPiDO -functionality. - -The example can be executed by running two files: - - cupido_example_step1.m - Matlab script to simulate data and generate example NetCDF - files with leveling and GPS data. The user can change the - simulation settings in the file. - - cupido_example_step2.py - Python2 script to convert the data in the NetCDF files to - an optimal set of double-differences. The user can change - the selection settings (area of interest, period of interest, - technique of interest, etc.) - -(c) Delft University of Technology, 2017 - +CUPiDO - Connecting Undifferenced Points in Deformation Observations + +This directory contains an example setup to demonstrate the CUPiDO functionality. + +The example can be executed by running three files: + + cupido_example_step1.m + Matlab script to simulate data and generate example NetCDF files with leveling and GPS data. The user can change the + simulation settings in the file. + + cupido_example_step2.py + Python2 script to convert the data in the NetCDF files to + an optimal set of double-differences. The user can change + the selection settings (area of interest, period of interest, technique of interest, etc.) + + cupido_example_step3_example_modeling.m + Matlab script to apply the inversion of two parameters of the Mogi model (i.e. depth and rate of the volume change), by direct search of a user-defined search space. The user can change the simulation and search space settings. + + +(c) Delft University of Technology, 2017 + diff --git a/example/cupido_example_step3_example_modeling.m b/example/cupido_example_step3_example_modeling.m new file mode 100644 index 0000000..0f1a1b7 --- /dev/null +++ b/example/cupido_example_step3_example_modeling.m @@ -0,0 +1,184 @@ +%CUPIDO_EXAMPLE_STEP3_MODELING Script to apply the modeling (inversion) of +%Mogi model parameters based on the results of steps 1 and 2 +% file. +% +% This script apply the inversion of two parameters of the Mogi model +% (i.e. depth and rate of the volume change ), by direct search of a +% user-defined search space. The objective function of the inversion is +% the L2-norm of the misfit between DD observations and DD predicted by +% forward models. The user can change the simulation and search space +% settings in the settings section below. +% +% (c) Sami Samiei Esfahany, Delft University of Technology, 2017. + +% Created: 18 July 2017 by Sami Samiei Esfahany + +clc +clear all +close all + + +data=load('cupido_example_result1.mat'); +%data=load('cupido_example_result2.mat') +%data=load('cupido_example_result3.mat') + + +%% DD observation vector (y) +y = data.DD_OBS'; % in [m] +NumDD=length(y); + +%% extract the coordinates of from/to benchmarks for each DD observation + +% give coordinate to the SYN_BM +[temp,ind_SYNBM]=ismember(cellstr('SYN_BM'),data.BENCHMARKS(:,1)); +if ind_SYNBM + data.BENCHMARKS(ind_SYNBM,2)=num2cell(single(10000000000)); % x coordinate of the SYN_BM (given outside of the AOI) + data.BENCHMARKS(ind_SYNBM,3)=num2cell(single(10000000000)); % y coordinate of the SYN_BM (given outside of the AOI) +end + +[temp,BMind_from] = ismember(data.DD_TABLE(:,1),data.BENCHMARKS(:,1)); +DDx1=cell2mat(data.BENCHMARKS(BMind_from,2)); +DDy1=cell2mat(data.BENCHMARKS(BMind_from,3)); + +[temp,BMind_to] = ismember(data.DD_TABLE(:,2),data.BENCHMARKS(:,1)); +DDx2=cell2mat(data.BENCHMARKS(BMind_to,2)); +DDy2=cell2mat(data.BENCHMARKS(BMind_to,3)); + + +%% Extract covariance matrix of DD observations +Qdd=data.DD_COV_MX; +invQdd =inv(Qdd); +%% extract the date of from/to epochs for each DD observation +[temp,EPind_from] = ismember(data.DD_TABLE(:,3),data.DATES(:,1)); +DDt1=datenum(cell2mat(data.DATES(EPind_from,2)),'yyyymmdd'); + +[temp,EPind_to] = ismember(data.DD_TABLE(:,4),data.DATES(:,1)); +DDt2=datenum(cell2mat(data.DATES(EPind_to,2)),'yyyymmdd'); + +%% extract sensitivity vector for each DD +DDsens=double(cell2mat(data.DD_TABLE(:,5:7))); + +%% extract technique of each DD +DDtech=data.DD_TABLE(:,8); + + +%% define the search space for the Mogi model +% only two parameters are estimated (depth and the rate of volume change) +% the coordinates of the center of the poit source are assumed to be known +% a-priori (deterministic) + +x0=10000; % y coordinate of the point source [m] - fixed (deteministic) +y0=10000; % x coordinate of the point source [m] - fixed (deteministic) + +depth_range=[2000 8000]; % range of the search space for the depth [m] +Vrate_range=[-800000 -300000]; % range of the search space for the rate of the volume chage [m^3/year] + +depth_int=100; % interval of the search space for the depth [m] +Vrate_int=5000; % interval of the search space for the rate of the volume chage [m^3/year] + + +d=[depth_range(1):depth_int:depth_range(2)]'; +v=[Vrate_range(1):Vrate_int:Vrate_range(2)]'; +[dd,vv]=meshgrid(d,v); % serach space + +ddall=dd(:); +vvall=vv(:); +NumSearch=length(ddall); + +%% loop on forward modeling +% i.e.: for each elemnt of the search space, we evaluate the Mogi model and +% construct the corresponding DD observations. Then we compute the misfit (residuals) +% between the model and DD observations and compute the L2-norm of the +% residuals. The model with minimum L2-norm is selected as the solution. +% Note that the L2-norm is computed in the metric of the inverse of the +% DD covariance matrix. + +L2norm=NaN(NumSearch,1); + +for i=1:NumSearch + + Depth=ddall(i); + Volume_rate=vvall(i); + + y_model_UP= (((DDt2-DDt1)/365.25)*Volume_rate*(1-0.27)/pi * Depth) .* ... + (((DDx2-x0).^2 +(DDy2-y0).^2+Depth^2).^(-3/2) - ... + ((DDx1-x0).^2 +(DDy1-y0).^2+Depth^2).^(-3/2)); % DD Mogi forward model (UP) + + y_model_E= (((DDt2-DDt1)/365.25)*Volume_rate*(1-0.27)/pi) .* ... + ( (DDx2-x0).* ((DDx2-x0).^2 +(DDy2-y0).^2+Depth^2).^(-3/2) - ... + (DDx1-x0).* ((DDx1-x0).^2 +(DDy1-y0).^2+Depth^2).^(-3/2) ); %DD Mogi forward model (East) + + y_model_N= (((DDt2-DDt1)/365.25)*Volume_rate*(1-0.27)/pi) .* ... + ( (DDy2-y0).* ((DDx2-x0).^2 +(DDy2-y0).^2+Depth^2).^(-3/2) - ... + (DDy1-y0).* ((DDx1-x0).^2 +(DDy1-y0).^2+Depth^2).^(-3/2) ); %DD Mogi forward model (North) + + y_model=sum(DDsens .* [y_model_N y_model_E y_model_UP],2); % construct the final DD model + e=(y_model-y); %residuals or midfit + L2norm(i)=e'*invQdd*e; % L2-norm of the residuals (in the metric of inverse of the DD covarianxce matrix) + + progress=100*i/NumSearch +end + + +%% find the mumimum L2-norm solution +L2norm_grid=reshape(L2norm,size(dd)); +[minL2,ind_solution]=min(L2norm_grid(:)); + +%% Visualization + +% the true values for visualization +Depth_True=6000; +Volume_rate_True=-550000; + +%% search space +figure +imagesc(d/1000,v,log10(L2norm_grid)) +hold on +plot(Depth_True/1000,Volume_rate_True,'ko','markersize',15,'MarkerFaceColor','w') +plot(dd(ind_solution)/1000,vv(ind_solution),'g.','markersize',30,'MarkerFaceColor','g') +temp = get(colorbar,'YTick'); +colorbar('Ytick',temp,'YTickLabel',10.^temp); +xlabel('Depth [km]') +ylabel('volume change rate [m^3/year]') +title('L2-norm over the 2D search-space') +legend('True value','L2-norm solution') + +%% Covariance matrix of DD observations +figure +imagesc(1e6*Qdd) +title('DD covariance matrix') +k = colorbar; +set(get(k,'title'),'string','mm^2'); + +%% DD model vs. DD observations + +% compute DDs predicted by the solution model +Depth=dd(ind_solution); +Volume_rate=vv(ind_solution); + +y_model_UP= (((DDt2-DDt1)/365.25)*Volume_rate*(1-0.27)/pi * Depth) .* ... + (((DDx2-x0).^2 +(DDy2-y0).^2+Depth^2).^(-3/2) - ... + ((DDx1-x0).^2 +(DDy1-y0).^2+Depth^2).^(-3/2)); % DD Mogi forward model (UP) + +y_model_E= (((DDt2-DDt1)/365.25)*Volume_rate*(1-0.27)/pi) .* ... + ( (DDx2-x0).* ((DDx2-x0).^2 +(DDy2-y0).^2+Depth^2).^(-3/2) - ... + (DDx1-x0).* ((DDx1-x0).^2 +(DDy1-y0).^2+Depth^2).^(-3/2) ); %DD Mogi forward model (East) + +y_model_N= (((DDt2-DDt1)/365.25)*Volume_rate*(1-0.27)/pi) .* ... + ( (DDy2-y0).* ((DDx2-x0).^2 +(DDy2-y0).^2+Depth^2).^(-3/2) - ... + (DDy1-y0).* ((DDx1-x0).^2 +(DDy1-y0).^2+Depth^2).^(-3/2) ); %DD Mogi forward model (North) + +y_model_solution=sum(DDsens .* [y_model_N y_model_E y_model_UP],2); % construct the final DD model +e=(y_model_solution-y); %residuals or midfit +%% +figure +plot(y_model_solution*1000,y*1000,'r.') +xlabel('DD predicted by the L2-norm solution [mm]') +ylabel('DD observations [mm]') +hold on +plot(abs(max([y_model_solution*1000;y*1000]))*[-1 1],abs(max([y_model_solution*1000;y*1000]))*[-1 1],'k--') +grid on +title('DD model vs. DD observations') +axis equal + +