Matlab Users
Contents
Layered elastic Models
Email from Tom Fournier (June 17, 2008)
Jeff recently got some layered half-space codes for an isotropic point source and a planer dislocation that run in matlab. I have been playing around with them, so if you are interested here are some details that might make your life easier. The functions are very easy to use, so really hats off to Kaj Johnson who developed them.
There are 4 m-files (that I assume are in /user/local/matlab, if not I'll put them there)
LayeredGreens.m
MogiLayers.m
get_scaleN.m
get_scaleN_mogi.m
The real work horses are LayeredGreens and MogiLayers. The inputs for these two functions are the same except for the model parameter vector. The functions require as inputs; a model vector, station locations, depth of layers, relative Lame constants for each layer, and an optional scale ( default scaleN=1). The output of these functions is a familiar 3xN matrix of site displacments.
I wrote a simple script that calculates the Lame constants from a velocity profile, vel2elastic, (now in usr/local/matlab). This takes as inputs the p-wave velocity, density, poisson's ratio, and optionally the p-wave/s-wave ratio. The density can be constant or vary by layer. Poisson's ratio is assumed to be constant unless you give the p/s ratio, then poisson's ratio is ignored and the elastic parameters are calculated using Vs. The outputs are Lame's parameters; Lame1 (lambda) and shear modulus.
ScaleN controls how many terms are used in the Hankel transform, which you can determine using get_scaleN and get_scaleN_mogi. These two functions simply compare the layered half-space solution to the homogeneous half-space solution. The output are the deformations from each case and a figure ploting the displacements for comparison. These functions call mogi.m and disloc3d.m, so make sure you have a path to a working version of these defined. The inputs are; a model vector, a maximim station distance, and scaleN. The maximum distance is so that you can compare results at the distance of your farthest station (or the farthest you want the solution to be accurate). ScaleN controls the accuracy of the solution but also the time it takes to calculate it, so depending on your problem, you want to pick scaleN large enough to give good results, but also small enough that you don't have to wait forever to get some answers. Ideally an optimization function should be written that will return a scale value optimized for time with a given tolerance between the homogeneous and layered half-spaces solutions. This type of function wouldn't be hard to write, nevertheless it does not exist in our lab.
A couple of caveats about get_scaleN functions. For a given scaleN, the accuracy of the layered model is dependent on depth, dip and other parameters that affect the depth to any portion of the fault. Keep this in mind when choosing a scale value.
The other warning I should bring up is in regard to the LayeredGreens.m function. This function will crash in cases where the layer depths are too deep. I simply changed the length scales from 1km to 10km, it worked sometimes. If you are interested in this problem, particularly if you want to try and fix it, let me know and I can tell you more about it.
Example: Creating a layered model
This example uses the velocity model at Veniaminof volcano [J.J. Sanchez, 2005, UAF thesis].
The information from this table is imported into MATLAB and converted into an elastic model using VEL2ELASTIC.m.
The variable "data" is the same as the table.
VEL2ELASTIC takes as inputs the p-wave velocity (km/s), density in each layer(kg/m3), and Poisson's ratio. Alternatively VEL2ELASTIC can also use the Vp/Vs (p-wave/s-wave velocity ratio), in which case it ignores the value of Poisson's ratio. The output are the Lame parameters for each layer.
The scripts that calculates deformation in a layered half-space, LayeredGreens.m (Okada dislocation), and MogiLayers (isotropic point source), use the Lame constants and the depth to the bottom of each layer. The bottom layer extends to infinity, so the depth vector has a length n-1, while the Lame parameter vectors have a length of n, for a model with n layers. **The layered model scripts require the length scale to be km.**
First make a vector with the depths to the bottom of each layer. Create a model vector, calculate the elastic parameters, and choose a scale term. **Note the scaleN is not required the default is 1.** The elastic parameter inputs need to be relative to the first layer, so divide by the first value in each vector. Notice that the program crashes, even though all of the inputs seem to be valid. You can get around this problem by adjusting the length scale, instead of using 1km as the length unit, use 10km. Now everything works just fine! Notice that the fault slip/opening does not have to be rescaled. Just remember that the units of displacement are the same as the units of slip/opening, whatever unit you choose.
Choosing a proper scaleN value is important, particularly if you are going to make many function calls to the layer model scripts. Use get_scaleN.m and get_scaleN_mogi.m to find the best scale value. The scaleN scripts require a model vector, a maximum distance, and a scale value.
>>m = [3 3 5 -45]; % For the Okada model only specify the length, width depth and dip >>rmax = 50; >>scaleN = 0.3; >> [Uprop,Uokada]=get_scaleN(m,rmax,scaleN);
The figure is created automatically, it is useful for double checking that the geometry of your model is what you think it is, and to be sure that the layered model solution matches the homogeneous solution. The best approach is to make scaleN as small as possible while still matching the analytical solution to within some tolerance, say 0.1mm or 0.01mm. Remember that the accuracy of the result, with a given scaleN value, may depend on the depth and dip of the dislocation, so if you are inverting for source geometry you should test several model geometries to ensure that you are picking a scale value that gives accurate results for all conditions.
Simulated Annealing
There are a handful of Matlab functions that are designed for easy implementation of the simulated annealing algorithm (Cervelli et al., 2001). The functions are located in /usr/local/matlab/TomSimAnneal/
The main workhorses are simanneal_gs.m and cool_simanneal.m while the other functions in the directory are useful for auxiliary tasks or play minor roles in process.
simanneal_gs.m
This is the actual simulated annealing algorithm, it is called _gs because of the geometric sampling method it uses (this isn't really important, but I thought you might ask).
The help output (type simanneal_gs at the command line) looks like:
% [allp MSE bestMOD lowMSE]=SIMANNEAL_GS(m,x,data,cov,CoSch,model) % % Simulated Annealing Algorithm model parameter inversion. Uses a geometric % sampleing of the model space to improve performance. % % %INPUT % m = model bondaries(4nx2),4 parameters for n isotropic sources, % min-max for each (east, north, depth, volume) % x = data locations (3nx1) % data = displacment/velocity (3nx1) % cov = data covariance (3nx3n) % CoSch = cooling scheduele (from cool_simanneal) % model = @function (w/ inputs: m,x) % -m = model parameters (vector) % -x = model inputs (site locations)(vector) % %OUPUT % PARAM = final model parameters(10xn) % MSE = vector of MSE values. % bestMOD = best fit model parameters (10xn) % lowMSE = lowest value from MSE, corresponds to model, bestMOD. % % Optimizations come from Cerevelli et al. 2001, JGR. % Extensive references regarding Simulated Annealing can be found in the % above paper. % % Tom Fournier; revised from sammeal_gs_sp Jan 2008.
This function inverts a function that you define in model. The only requirement is that model accepts only two input vectors; m the model parameters, and x additional model inputs (for GPS this is often site locations).
bootstrap.m layermogiform.m nlinfit_cov.m
layerokadaform.m scatter_sample.m
