Although its getting a little old, the primary documentation is contained in the TeX file
(TeX is an ascii-based typesetting code, so this file is actually readable after a fashion without compiling it with TeX, although free, public domain versions of TeX are available for all platforms). Here's a HTML version of the document:
If you are in the business of MT inversion, you are encouraged to consider the `Geotools' by Geotools Corporation, which embodies Occam inside a nice interpretation and data management package. Geotools sponsored a lot of the version 2.0 re-write which makes up the current release and so are a worthy crowd.
You are free to pass Occam on to third parties AS LONG AS YOU HAVE NOT ALTERED THE CODE.
If you, or the third party, let me know I will add them to the distribution list for bug reports, upgrades, etc. If there really is some reason you want to alter the code, let me know and I will consider implementing your improvements for all.
Some of the files that you will encounter are:
Document.TeX: This document containing release notes, occam documentation and mt2D documentation
occam.f: Contains the inversion routines needed for any 1D, 2D or 3D inversion.
occamdim.inc: Array dimensions for occam.f
occam.inc: Common blocks for occam.f
imp.inc: Included in every subroutine for general compiler directives.
mt2D.f: Contains the routines needed by occam.f to perform 2D MT inversions.
mt2Ddim.inc: The array dimensions for mt2D.f
mt2D.inc: Common blocks for mt2D.f
mt2DExample2: An example based on synthetic data generated over a conductive (10 Ohmm) prism in a 100 Ohmm background and 5\% noise in the data. Break this file up into MESH, INMODEL, DATA.TEST2, and startup.
MakeModel2Dmt.f: Poor man's Geotools. Makes the startup files needed to run occam2Dmt, given a data file and answers to some simple questions.
FindResponse2Dmt: finds the response of a model in an ITERxx file and also outputs the model as a set of postscript commands to draw coloured blocks.
SynMod2Dmt: Generates a set of synthetic data with noise from a startup or ITERxx file.
mt1D.f: Contains the routines needed by occam.f to perform 1D MT inversions. These routines were cobbled together from existing routines in the old OCCAM1 release and are not very neat. Any problems you find in these routines you should be able to fix yourself.
mt1D.inc: Common blocks for mt1d.f
MakeModel1Dmt.f: Program to construct input files for simple 1D MT inversion.
mt1DExample1: Example based on the COPROD1 data set.
PlusD.f: Bob Parker's 1D MT D+ code, included with permission.
The Occam, mt2D, and mt1D documents below should tell you how to get started.
Document Revision 1.01, Occam revision 2.01
Introduction.
The OCCAM code has undergone major revision as part of the complex job of performing 2-dimensional (2D) magnetotelluric (MT) inversions with a variety of data, extra constraints, extra or relaxed penalties, and fixed or additional parameters. The fundamental engine of the inversion is unchanged, with the exception of an added penalty against a preferred model, or prejudice. It has, however, been extensively reorganized in order to decouple it completely from the forward problem being inverted, and to take some advantage of dynamic memory allocation or memory overlay. This reorganization improves one's ability to use OCCAM for problems other than MT.
My approach to adding more complications to the 2D MT inversion, done at the request and under the sponsorship of WSE Associates, has been to completely purge OCCAM itself of the nasty problems I was faced with. The only cost from a user's point of view of this rather cowardly action is that the generation of the penalty matrix is now part of the responsibility of the forward problem. While the traditional and simple 1D penalty matrix is easily generated using the subroutines provided, by removing this job from OCCAM we gain substantial flexibility in tailoring problems to our needs without having to re-write the inversion routines. It turns out that the penalty matrix is the only difference between 1D and 2D problems from OCCAM's point of view. Occam's world is composed of simple 1D vectors of data constraints and model parameters, linked together by the 2D Jacobean matrix. The penalty matrix simply takes pairs of model parameters and adds a weighted difference to the penalty measure. Traditionally, this job has been simply to penalize adjacent pairs of 1D model layers, or vertically and horizontally adjacent pairs of 2D model bricks. However, any pair of model parameters can be penalized, as often and as much as desired or not at all. Each additional penalty simply becomes an additional row of the penalty matrix, and since OCCAM only wants to use a matrix which is the number of parameters (NP) square, it does not matter how big the penalty matrix gets.
The notable change to the inversion engine is that I have taken the algebra that Catherine deGroot-Hedlin used to penalize the free parameters into peripheral fixed structure and generalized it to provide a penalty into a preferred model consisting of a complete set of model parameters. OCCAM always does this now, since the algebra requires modification deep in the heart of the beast, but usually the weights for this penalty will be set to zero. Of course, strictly speaking it is the user's responsibility to set these weights.
As the forward problems got slower and more complicated, the desire to sit and talk to the inversion program diminished, and it became important that restarting of the program from the middle of an inversion was made easy. Also, as the inversion became truly model independent, a standard set of inversion parameters could be supplied. Consequently, I have written a calling routine which reads in all the relavent information from a file, called 'startup'. 'startup' points to a data file and a model file, which may point to other files in turn.
Of course, all these files have to be constructed, which will usually be done by the very computer code that was built into the 1D and 2D versions, so again we have removed this problem-dependent nastiness and preserved the sanctity of inversion.
Portability
The routines are written in FORTRAN 77. An internal list-directed read in version 2.0 has been removed so that the program will run on the Unix Cray as is. Most of the subroutines are pretty bullet-proof by now. I have converted everything to single precision, to keep memory requirements down for large problems. Your forward routines are welcome to use double precision, as do the 1D MT routines I supply, but make sure you pass out single precision version of the final response. No penalty associated with going to single precision has become apparent. All declarations are expressed as REAL or INTEGER, and many compilers will let you switch the default precision using compile-line options, so if you have any doubts or paranoia and a lot of memory, try compiling for extended precision.
The user will undoubtedly want to use common blocks to convey information from INPUTM() to FORMOD() and OBJMAT() (see below). Since these commons are not be declared in the calling program, you need to make sure that your compiler uses static commons or can be set to do so. Static commons will speed up execution a lot, anyway. Local variables are assumed to be dynamic, and anything that needs to be saved is explicitly SAVEd.
I have tried to cover the use of the code in these notes and in the notes that accompany the specific inversion packages. If you have any problems or questions, look at the code before you call me. It is complicated but not obscure, and has a great deal of self-documentation.
Running OCCAM.
a) First check the maximum dimensions of the problem, which can vary tremendously, in the include file `occamdim.inc'. Since these are maximum dimensions only, they do not have to be changed unless your particular problem will exceed them. Normally there will be additional dimensions associated with fixed model parameters (layer thicknesses, brick sizes, etc.), held in user's common blocks for use by the forward routines, and attempts to increase the size of the model should check these as well (e.g. for mt2D check the dimensions in `mt2Ddim.inc'). There might also be dimensions associated with the size of the data array (e.g. maximum number of stations allowed in 2D inversions) in user's commons. However, all that Occam needs is
For simple 1D, NNP = NPP, while for simple 2D NNP a little less than 2*NPP. Bear in mind that any embellishment of the simple penalty matrices, such as deleting or adding penalties, will vary these sizes.
b) Six subroutines have to be supplied to compute the forward response of models and derivatives, input the data and starting model, compute the penalty matrix, and provide a date/time stamp for output files. Before I describe these routines, I should mention the three common blocks used by Occam, which are contained in `occam.inc':
common /data/ nd, d(ndd), sd(ndd), dp(ndd,4), ndp
/data/ contains the constraint data, data errors and some information about the data. The user-supplied subroutine INPUTD() is responsible for setting the contents of this common. The user's forward modelling routines FORMOD() and FORDIV() will be passed nd and dp() through the argument list.
common /model/ np, pm(npp), dm(ndd)
- np [integer] is the number of free model parameters
- pm() [real] is the array of model parameters
- dm() [real] is an array used to carry model response (data domain) around
/model/ contains the current model parameters for this iteration and its response. The user does not actually have to set any of this, as Occam will fill these arrays using parameters read from the startup file.
common /result/ ptp(npp,npp), wjtwj(npp,npp), wjtwd(npp), pwk1(npp), pwk2(npp), pwk3(npp), frac, nfor, iout, idebug, prewts(npp),premod(npp)
/results/ contains stuff to be carried into subroutines of OCCAM. Most of the stuff is of no concern to the user, but if a preferred model is to be used the premod() and prewtd() vectors need to be set to the preferred model vector and penalty weights, respectively. The penalty weights are initialized to zero by OCCAM, so if the user does nothing then nothing will happen.
And so on to the routines which you have to supply:
FORMOD(NP,ND,PM,DP,DM): This routine produces the model response (synthetic data) given a vector of model parameters. The user is likely to want to make common blocks available to FORMOD() which contain fixed model parameters (layer thicknesses, block geometry, etc.) which have been set up by INPUTM(). This is fine, but don't attempt to make /model/ available, as the inversion relies on being able to pass various attempts to find a good model in to FORMOD() via PM(), while maintaining the current model in /model/ PM(). Similarly, since DP() and ND are passed in through the argument list, don't try to include /data/ either.
On input:
On output:
- DM(ND) (real vector) is a vector containing the model response, in the order dictated by DP() and so in one to one correspondence with the real data in D().
FORDIV(NP,ND,PM,DP,DM,WJ): This routine is the same as FORMOD except that the Jacobean or sensitivity matrix, WJ(), is computed as well as the forward response. WJ(i,j) is the partial derivative of the {\it i\/}th datum with respect to the {\it j\/}th model parameter.
OBJMAT(IRUF, DEL, NPEN): Is a routine which generates the penalty matrix. This routine is called twice every iteration, once to set up PTP()=$\partial^T\partial$, and once to compute the model roughness, so don't destroy anything you need on the first call. Again, the user will probably want to use INPUTM() to set up common blocks containing information about exceptions to the default penalty matrix (such as jumps between layers or blocks) which will be passed to OBJMAT().
On input:
On output:
DATIME(DATETM): Is a routine which returns an 80 character string (DATETM) containing the date and time the subroutine was called. This routine will be system dependent, and can always be written to return a string constant declaring that the date and time are not available. This is what the distribution code does.
INPUTD(FILENAME): Is a routine which takes the name of a file (FILENAME) and reads from this file the data required for the inversion. It is the responsibility of this routine to set ND, D(), SD() and DP() in /data/.
INPUTM(FILENAME): Is a routine which takes the name of a file (FILENAME) and reads from this file all the information about the model that will be required by the forward routines FORMOD() and FORDIV(), as well as information required by OBJMAT(). This information must be stored in the user's own common blocks shared by FORMOD, FORDIV, OBJMAT, and INPUTM.
c) The startup file needs to be created. This file, called `startup', contains the starting model vector and other inversion parameters. At each iteration a version of this file will be generated, called ITERxx where xx is the last iteration number. The format of the ITERxx files is identical to the startup file, so any of these can be renamed `startup' in order to continue the inversion where it left off (if you do this do not alter anything except the file's name, and make sure that the user's data and model files are still available). So that the original startup parameters are not lost they are copied immediately into ITER00 before the first iteration.
Most of the lines in the startup file begin with an 18 character label stating the nature of that line's contents. This is purely for readability of the file, as these labels are skipped by OCCAM. There are two consequences of this, however. Some character strings must start exactly in column 19 before they are read as intended, and OCCAM uses the equivalent of an internal free-format read from the string starting in column 19 for all integers and reals that are required.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only OCCAMITER\_1.0 is supported. Beware of leading blanks after the 18 character label.
Field 2: 18 character `DESCRIPTION' label followed by 80 character description. This can be any information the user wishes to attach to the file. Ignored by OCCAM.
Field 3: 18 character `MODEL FILE' label followed by 80 character model file name, which will be passed into INTPUM(filename).
Field 4: 18 character `DATA FILE' label followed by 80 character data file name, which will be passed into INTPUD(filename).
Field 5: 18 character `DATE/TIME' label followed by 80 character date stamp. This is ignored by OCCAM and is for the user's benefit only. Anything can be written here in the startup file, and the ITERxx files will have a valid date and time stamp or a fixed string, depending on what system-dependent routines are available.
Field 6: 18 character `MAX ITER' label followed by the maximum additional iterations to be performed this run. That is, inversion stops after iteration (ITERATION + MAX ITER) is performed.
Field 7: 18 character `REQ TOL' label followed by the desired root-mean-square misfit that the smooth model's response is to produce.
Field 8: 18 character `IRUF' label followed by the type of roughening operation required (0,1,2). This integer is passed to OBJMAT(). Note that your OBJMAT() routine does not have to support all of (0,1,2), but there is no point in supplying unsupported options here. OCCAM does not otherwise use IRUF.
Field 9: 18 character `DEBUG' label followed by the debug level. 0 is silent, 1 is just as fast but verbose, 2 produces a blanket search of lagrange multiplier which is slow but more robust than the normal operation.
Field 10: 18 character `ITERATION' label followed by the number of the current iteration. This should be set to 0 in the original startup file and will then be set by OCCAM.
Field 11: 18 character `PMU' label followed by the log10 of the Lagrange multiplier used for the current iteration. This should be set to any readable number on the zeroth iteration. It will be set by OCCAM on subsequent iterations and should be preserved at this value when restarting.
Field 12: 18 character `RLAST' label followed by the roughness of the current model. This should be set to any readable number on the zeroth iteration. It will be set by OCCAM on subsequent iterations and should be preserved at this value when restarting.
Field 13: 18 character `TLAST' label followed by the misfit of the current model. This should be set to any readable number on the zeroth iteration. It will be set by OCCAM on subsequent iterations and should be preserved at this value when restarting.
Field 14: 18 character `IFFTOL' label followed by an integer that indicates whether the current model is feasible or not (1=yes, 0=no). This should be set to any readable number on the zeroth iteration. It will be set by OCCAM on subsequent iterations and should be preserved at this value when restarting.
Field 15: 18 character `NO. PARMS' label followed by the number of model parameters in the inversion, NP.
Field 16: NP entries (free format) of model parameters. A starting model must be provided, usually perfectly smooth but having some realistic value. Normally this model image would not be altered during a re-start. There might be circumstances in which one might want to manually change this model, but since PMU, RLAST, TLAST and IFFTOL are all tied to the model appearing in this field, if it is changed it would be best to set ITERATION to zero to start the inversion from scratch.
d) Compile and link occam.f with the user-supplied routines DATIME, FORMOD, FORDIV, OBJMAT, INPUTM, INPUTD, and make a startup file available. Provide the data file and model files required by your input routines INPUTM and INPUTD, and run the program.
If You Want To Know More, or, If Anything Goes Wrong
Data/Model Size. The distribution code is set up to handle up to 200 data and 20 model parameters. Make sure that these parameters are set large enough for your problem (see` occamdim.inc'), and that any arrays which you are using for the forward model and object matrix computations are also large enough.
There is a little bit of error checking inside OCCAM, and, for example, if you try to read too many parameters from the model file you will be told of your folly. When your write the subroutine OBJMAT(), you should include OCCAMDIM.INC to get the dimensions you need, and it is a good idea to check the size of the penalty matrix against NNP, line by line, so that you can stop before you go too far.
Forward routines. Specifications for the forward routines FORMOD() and FORDIV() are given above and in the code. Like all linearised methods, errors in the derivative matrix may lead to problems. It is best if the derivatives are computed from analytical expressions (remember that filter or fast Hankel transforms generally allow the kernel to be differentiated and then filtered to give the derivative; see the appendix of Constable et al., 1987). Always check your routine against a finite difference calculation before using it for inversion. Don't forget that if data fitting is done in log domain and/or layer resistivities (or conductivities) are expressed as logs then the derivatives must be appropriately scaled.
Here's a bit of code that will use central differences for FORDIV(): fordiw.f
Lagrange multiplier and expected fit. Both the strength and weakness of OCCAM lies in choosing the Lagrange multiplier. It relies on univariate minimization and root finding schemes that can get confused by local minima, and their performance will really depend on how non-linear your forward model is. Running with debug level 2 will scan 8 decades of Lagrange multiplier looking for a global minimum and will be helpful if you have real problems, but this will require two or three times the number of forward computations per iteration that normal operation does.
The choice of target R.M.S. misfit is important. Attempting to over fit the data by requesting too small a misfit will fail. Minimum second derivative modelling seems more sensitive in this respect than minimum first derivative modelling. Try and get a feel for what represents a reasonable fit to your data before choosing a target misfit; after all, obtaining an expected, rather than best, fit to the data is the whole philosophy of the method. For M.T. and resistivity, Bob Parker's D+ and bilayer algorithms will provide the best possible fits obtainable with a 1D earth. Don't expect to do as well as these with a smooth model! Simulations suggest that fitting about 25\% worse than the best fitting model is reasonable. If an existence algorithm such as D+ is not available, OCCAM itself will get quite close to the best fit in many circumstances. Run OCCAM towards a low misfit, and then rerun it to obtain a smooth model for a misfit about 1.25 times the best you can get. (I have included, with Bob's permission, the source code for D+ in the 1D MT package). plusd.f
Adequate model space. If the program cannot attain what you consider to be a reasonable fit to the data, make sure that the model space being used has enough freedom to model the data. For 1D EM inversion there should be at least 4 samples of the model (or layers, if a layered structure is used) per decade. More is better, especially if you are trying to find the best fit possible on an initial investigation of the data. Examine the output model to make sure that there is no structure at the edges (in the deepest and shallowest parts for 1D). Going about a decade deeper and shallower than the limits of resolution of the data with the model is quite reasonable, even if the associated depths look ridiculous.
Remember, there is a penalty on model structure and so you are not going to get into trouble by `over-parameterizing' things. More is better as far as OCCAM is concerned (computer time and memory are your problems, not Occam's).
Convergence. Firstly, I want to say that OCCAM should converge (step size goes to zero), and that you should let it do so. It might appear that once one has achieved the desired misfit further use of computer time is wasteful, or at least delays you getting your hands on the model, but this is not so. It takes several iterations after the desired misfit has been attained to remove unnecessary structure from the data, and you do not have the truly smoothest model until OCCAM has converged. If you do not want the smoothest model, there are probably algorithms out there that will find you a feasible model much faster than OCCAM will. If you do want smooth models, then let there be convergence.
A step-reduction algorithm attempts to force convergence when things go bad. Bad in this case is failing to find a minimum misfit that is smaller than the last one, or, after the required misfit is attained, a roughness that is smaller than the last one. The step-reduction algorithm usually gets things by but is not very good.
Convergence to perfectly smooth models (half spaces for 1st derivative and ramps for 2nd derivative smoothing) can be troublesome, usually because there is more than one solution. OCCAM1 derives a great deal of stability from the uniqueness of its solution, and may play up if there are two equally smooth models which both fit the data equally well, as there will usually be if one perfectly smooth model fits the data adequately. The program keeps an eye on the roughness, and if it gets so small it thinks the model is perfectly smooth and the desired misfit is attained then it will stop and declare a finish.
For similar reasons, the main mode of failure in OCCAM occurs when the program ends up jumping between two models, both of which have similar misfit and roughness and which `point' to each other. I have not attempted to fix this yet.
If you have convergence problems and can afford to run the forward program a hundred times or so, setting DEBUG to 2 in the startup file will print a table of misfit versus Lagrange multiplier. This is very helpful in assessing the behaviour of the forward problem and the performance of the inversion. Furthermore, the code will use all this extra information to avoid local minima.
One final note on convergence. Bob Parker cannot work out why OCCAM works or prove that it will ever converge to anything.
Catherine deGroot-Hedlin pioneered the work on extending the regularized inversion to the 2D magnetotelluric (MT) problem. Her thesis contains all the work she did on this subject, but the relevant stuff is contained in 2 published papers and one more that is still in progress as of this date. Under sponsorship from WSE Associates, Catherine and I tidied up the code after her thesis work was finished. This finally resulted in rewriting most of the routines from scratch.
As well as a list of functions that WSE wanted implemented, I had a driving philosophical desire to remove every single aspect of the forward modelling routines from the inversion routines. This also has the practical advantage that the inversion routines can be thoroughly tested on the much faster 1D case and extension to other problems and 3D.
Those of you who have Catherine's original code have a functionally working system. However, that code could not be maintained in Catherine's absence, and it was also quite complicated to use. The new code can do quite a few things that the old code could not do (notably to include topography), but the main thing is that the file structure it uses allows more complicated tasks to be performed almost as easily as a vanilla inversion. However, for the examples that appear in Catherine's papers and thesis, it produces the same results as the original code.
Our version of the 2D magnetotelluric (MT) inversion relies on Phil Wannamaker's finite element (FE) forward code. This code, although slow, was the only publicly available code that includes calculation of the derivatives needed to construct the Jacobean matrix, necessary for linearized inversion. We acknowledge Phil's cooperation in this regard; this work would never have happened without his code.
Distribution and Portability
The 2D MT package consists simply of files `mt2D.f', `mt2Ddim.inc', `mt2D.inc'.
The routines are written in FORTRAN 77, and I explicitly take advantage of the FORTRAN 77 feature that states that thou shall not necessarily go through a do-loop. The common blocks used to convey information from INPUTM() to FORMOD() and OBJMAT() are not declared in the calling program, so you need to make sure that your compiler uses static commons or can be set to do so. Local variables are assumed to be dynamic, and anything that needs to be saved is explicitly SAVEd. The code has currently been run on a Macintosh, Spark 2 and a Unix Cray.
An overview of using occam2Dmt.
a) Set the maximum dimensions for your problem, if they are not already adequate. Dimensions associated with the forward code are located in `mt2Ddim.inc', those associated with the inversion in 'occamdim.inc'.
b) Compile and link occam.f and mt2D.f, making sure that all the include files are available.
c) Create the input files. This will normally be done by another computer program, such as WSE's Geotools or the simple program that I distribute with the academic version of the code. occam2Dmt does only minimal checking to make sure the various inputs are valid and consistent, relying on the code that generates the files to get it right.
There are 3 to 5 files required to run the occam2Dmt, in addition to the startup file required by OCCAMv2.0. These are the data file, model file, FE mesh file and the optional static shift and model prejudice files.
Setting the dimensions
The Occam distribution notes describe setting NPP, NDD and NNP in occamdim.inc. Basically, these numbers are the maximum number of model parameters (regularization bricks), maximum number of data (sites x frequencies x components), and maximum number of penalty terms. In the distribution, NNP is set to 2*NPP, which will accommodate vanilla 2D inversion.
The dimensions for the 2D MT problem are stored in the file `mt2Ddim.inc'. They are
IPY = Maximum number of nodes in the FE mesh in the Y direction (NODEYin the MESH file). Note that nodes = elements + 1.
IPZ = Maximum number of nodes in the FE mesh in the Z direction (NODEZ in the MESH file). When computing TE responses, IPZ must exceed NODEZ by at least 10.
IPI = Maximum number of regularization bricks in the inversion model (`NO. PARMS', or NP, in the startup file). The same as NPP in occamdim.inc.
NL = Maximum number of layers on the left or right side of the FE mesh, that is, the number of regularization brick in the vertical direction (`NUM LAYERS', or NLAYER, in the MODEL file).
NF = Maximum number of frequencies (`FREQUENCIES', or NFRE, in the DATA file)
NR = Maximum number of stations (`SITES', or NRC, in DATA file)
MCOL = Maximum number of columns including layering at sides, that is, the number of regularization blocks in the horizontal direction (the largest value of NRCOL(I) in MODEL file).
maxexp = maximum number of additional rows to regularization matrix (`NO. EXCEPTIONS', or NEXEP, in MODEL file).
The Data File:
This file is pointed to by the startup file, but points to no other files. In the description below, a `field' is usually one line, but not always so.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only `OCCAM2MTDATA\_1.0' is supported. Make sure there are no leading blanks after the 18 character label.
Field 2: 18 character `TITLE' label followed by 80 character title. This can be any information the user wishes to attach to the file. Ignored by occam2Dmt.
Field 3: 18 character `SITES' label followed by number of sites, NRC, in free format.
Field 4: NRC site names, one per line. These are ignored by occam2Dmt, and serve for user identification only.
Field 5: `OFFSETS' label. Skipped by occam2Dmt.
Field 6: NRC station offsets, free format. These are site coordinates in metres, positive towards the right of the model. The origin can be any arbitrary point; the position of the left edge of the model in the same coordinate systems appears as `BINDING OFFSET' in the model file. The site locations will be moved to the nearest FE node (under debug level 1 or greater, the nodal locations of the sites will be printed at the start of the inversion).
Currently, the FE forward subroutine PW2DI() requires that station offsets be in order of increasing value, and that no two sites may be on the same FE node. It is expected that this restriction will be lifted in the future. It is probably not a good idea to have a station on the edge of a regularization brick, but rather have them place on FE nodes within a brick. It is not clear to me how important a consideration this is, since the surface structure is usually smooth, but you have been warned.
Field 7: 18 character `FREQUENCIES' label followed by the number of frequencies, NFRE.
Field 8: NFRE frequencies, free-format, in Hz.
Field 9: 18 character `NO. DATA BLOCKS' label followed by the number of data
blocks (free format), NDATA.
Field 10: `SITE, FREQUENCY, DATA TYPE, DATUM, ERROR' banner. Skipped by OCCAM2DMT.
Field 11: NDATA lines of data (free format). Each line is made up of the site number (index of the station offsets in field 6), frequency number (index of frequencies in field 7), data type (see below), datum, and standard error. The data types are as follows:
The Model File:
This is pointed to by the startup file. It points in turn to a mesh file containing a valid description of a PW2DI (Wannamaker) FE mesh. This mesh file provides the yardsticks for the description of the regularization mesh. See deGroot-Hedlin and Constable (Geophysics 55, 1613-1624) for the distinction between the regularization mesh and the FE mesh. The model file also points to an optional static shift file and an optional model prejudice file.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only `OCCAM2MTMOD\_1.0' is supported. Beware of leading blanks aftet the 18 character label.
Field 2: 18 character `TITLE' label followed by 80 character title. This can be any information the user wishes to attach to the file. Ignored by occam2Dmt
Field 3: 18 character `DESCRIPTION' label followed by 80 character model description. This can be any information the user wishes to attach to the file. Ignored by occam2Dmt
Field 4: 18 character `MESH FILE' label followed by 80 character mesh file name. occam2Dmt will attempt to open this file to read the FE mesh description.
Field 5: 18 character `MESH TYPE' label followed by mesh format. Currently only `PW2D' is supported.
Field 6: 18 character `STATICS FILE' label followed by 80 character statics file name. occam2Dmt will attempt to open this file to read the static shift information. However, `none' or `NONE' are valid keywords which will simply set all static shifts as fixed parameters of size zero (i.e. not implement static shift corrections or inversion).
Field 7: 18 character `PREJUDICE FILE' label followed by 80 character model prejudice file name. occam2Dmt will attempt to open this file to read a set of preferred model parameters and weights. However, `none' or `NONE' are valid keywords which will simply set the weights for the preferred model to zero (i.e. not implement a preferred model).
Field 8: 18 character `BINDING OFFSET' label followed by (free format) `binding offset', BOFFS. This is the position, in metres, of the right edge of the left terminating layer in the regularization (not FE) mesh. This number essentially sets the origin for the problem and is used to convert station locations, in metres, to station locations in node number. This allows the original survey coordinates of the MT stations to be preserved.
Field 9: 18 character `NUM LAYERS' label followed by (free format) number of layers in regularization mesh, NLAYER. This number includes (in 1D terminology) the terminating half-space. Layers are counted from the surface down.
Field 10: This is the description of the regularization mesh in terms of the FE mesh. It consists of NLAYER repeats (I=1,NLAYER) of two subfields:
Subfield 10.1: IRZ(I), NRCOL(I) (free format). IRZ(I) is the number of finite elements making up the vertical thickness of the Ith layer, and NRCOL(I) is the number of regularization bricks across the Ith layer.
Subfield 10.2: NRCOL(I) occurances of IRY(I,J), the number of finite elements making up the horizontal width of the Jth brick (counting left to right) in the Ith layer. Negative numbers will flag the brick as being composed of entirely fixed structure.
This description can be used to make the regularization bricks increase in width as one gets deeper, keeping the aspect ratio approximately unity or whatever is appropriate for the problem. Currently, the only restriction on the nature of the change in size of the regularization bricks from layer to layer is that at least one corner of each brick must be in common with a corner of the bricks in the layers above it, and similarly a corner must be shared with a corner below it. This is because the algorithm that searches for adjacent bricks during the construction of the penalty matrix looks for matching corners of adjacent bricks. This is not a severe restriction, but it is expected that it will be lifted in the future.
The entire FE mesh (see below) must be overlain by the regularization mesh, even if there is fixed structure in the model. Any regularization bricks that are entirely fixed structure will be dropped from the model parameterization.
Field 11: 18 character `NO. EXCEPTIONS' label followed by (free format) number of modifications required to the default penalty matrix, NEXEP. NEXEP=0 is allowed. The default penalty matrix penalizes all horizontally adjacent pairs of free regularization bricks and all vertically adjacent pairs of free regularization bricks. The penalties increase exponentially with depth (although currently this is enforced through the assumption that the layer thicknesses increase exponentially with depth). Modifications to this scheme can be entered here and fall into two types: changing the weight the default entries (to zero in order to relax the penalty at a known discontinuity or to a large number to lock two bricks together) or adding additional penalties between non-adjacent bricks. In either case, the scheme works by you providing the layer and column indices (in units of regularization bricks) of the two bricks to penalize followed by a penalty weight, which is a multiplier of the default penalty for the depth associated with the first brick:
Field 12: NEXEP free format entries of L1, C1, L2, C2 [all integer], PEN [real]. (L1,C1) are the layer and column number for the first brick of the pair (I and J in the notation of Subfield 10.2), similarly (L2,C2) for the second brick, and PEN is the penalty to apply between them.
The Mesh file:
This is a standard PW2DI model description. It is identical to Phil's description with only two exceptions: 1) an extra alphanumeric is allowed in addition to 0 (air), 1-10, A-Y (all fixed resistivities), and Z (seawater). The added character is `?' and indicates that a free parameter associated with the regularized inversion is to be substituted in this location. In this way fixed and invertible structure can be defined as part of one FE mesh. There is no restriction on where fixed structure can occur; each triangular FE sub-element is examined individually for fixed structure. In this way part of a regularization brick can be fixed structure and the rest becomes a parameter. Topography and bathymetry can be described on the upper surfaces of regularization bricks, since air and seawater are fixed structure.
2) IDX, the option type, has been extended to include 0, indicating compatibility with Occam. For those of you who do not have Phil's document, here is an abridged description of the mesh file. In any situation in which Phil's description duplicates parameters given in the other Occam2DMT data or model files (e.g. number of frequencies, data type), the other files take precedence and the value in the mesh file is ignored. However, if you keep these numbers sensible, then you will be able to run the mesh file with Phil's code, independently of Occam. This gives you more versatility and keeps me honest.
Field 1: 80 character title or description, ignored by Occam.
Field 2: 6 free-format integers:
Field 3: NRES (up to 35) resistivities associated with fixed structure in the model, not including air. This forms the lookup table for substitution of the element flags set in Field 8 (below), and must occur in the order of these flags (1-9, A-Y). Free format.
Field 4: NFRE frequencies at which model response is calculated. Ignored by Occam. Free format.
Field 5: NODEY-1 FE (column) widths, metres, free format.
Field 6: NODEZ-1 FE (layer) heights, metres, free format.
Field 7: NRC = Number of stations (ignored by Occam, can be 0), followed by NRC entries of nodal positions for each station (also ignored by Occam).
Field 8: Array of alphanumeric labels for each element of the FE mesh. Format is (170a1) (i.e. do not exceed 170 entries on each line). The labels go from 1-9 and A-Y, corresponding in order to the entries in Field 3 (above), and `?', which indicates a free parameter for the Occam inversion. Air is 0 and sea is Z.
Each finite elephant is divided into four triangular sub-elephants by two diagonal lines, giving a top, left, bottom and right triangle. Horizontal ordering (columns) of labels corresponds to horizontal elements, so normally each line would be NODEY-1 labels long, although the lines can be broken into two or more as long as the last horizontal element comes at the end of a line. The vertical (row) ordering is firstly by triangular element, in the order of top, left, bottom, right, and secondly by vertical element position. Thus the first four lines would normally be the four triangular sub-elements of the first row of elements, the next four lines the four sub-elements of the second row of elements, etc. There will normally be 4 x NODEZ-1 records in this field.
This is the end of the mesh file, although Phil defines some more parameters after field 8 for polygonal inversion. These parameters may be there; they will just be ignored by Occam.
Static shift file:
The optional static shift file provides for fixed or invertible static shifts. Fixed static shifts appear in this file and will be applied during every forward calculation. Invertible static shifts are defined as such in this file, added as parameters to the array PM after the model resistivity bricks have been taken care of, and will take on initial values from the parameter vector in STARTUP (so the user, or the user's model generation program, has to keep track of how many parameters are needed and of what type they are). The values given in the statics file are not necessarily ignored, however. If a data constraint is required to control the magnitude of the average free static shift, that constraint is taken from the sum of the values associated with the free static shifts in the statics file.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only `OCCAM2MTSHIFT\_1.0' is supported. Beware of leading blanks after the 18 character label.
Field 2: 18 character `DESCRIPTION' label followed by 80 character description. This can be any information the user wishes to attach to the file. Ignored by OCCAM2DMT.
Field 3: 18 character `DATA FILE' label followed by 80 character data file name. This is provided only so that the user can track the relationships between the data and static shifts, and is ignored by OCCAM2DMT, which takes the data file name from the STARTUP file.
Field 4: 18 character `CONSTRAINT TYPE' label followed by the free format integer ICONS. ICONS controls the type of data constraint used to control the baseline of any free static shifts (if all the static shifts are fixed, ICONS should be set to zero, although OCCAM2DMT will do this for you).
ICONS=0: No data constraint is used to control the static shift baseline. This option would be used when the model resistivity is controlled by penalizing into a model prejudice, or when a few stations with fixed static shifts are included. Also, free static shifts can be associated with preferred values and penalized directly.
ICONS=1: The sum of the combined free TE and TM static shifts (in log domain) are forced to sum to the same value as the sum of the values appearing in this file. Note that the sum is computed only for the free (invertible) shifts and does not include shifts included as fixed parameters. Remember that this data constraint contributes to a 2-norm misfit measure; be careful in setting the error term (below). Normally, the shifts will be required to sum to zero, so all free static shifts will be set to zero in this file.
ICONS=2: Same as ICONS=1 except the free TE and TM static shifts sum separately to their starting sums.
Field 5: 18 character `CONSTRAINT ERROR' label followed by the free format integer SSERR. SSERR is an error used with the data constraints type 1 and 2 to allow some flexibility in applying the constraint. Try 0.1 to start with (sum of free static shifts will come to within about 0.1 log unit of starting sum).
Field 6: 18 character `NO. SHIFT BLOCKS' label followed by the free format integer NSB. NSB predicts how many blocks of static shift data appear below.
Field 7: Descriptive banner ignored by OCCAM2DMT.
Field 8: NSB repeats of static shift information: integer site number, integer data type, real shift value, integer fixed/invertible flag.
The site number and data type are identical to those used in the data file. Only data types 1 (TE resistivity) and 5 (TM resistivity) are currently supported, although it is intended that type 0 will be added to indicate a shift which is the same for both TE and TM modes at each site. You do not have to represent all sites and data types, as static shifts for any unrepresented data will be set to fixed and zero. Shift is the log10 value of the multiplicative shift at each station/mode. That is, the additive distance on a log10 plot over which the curves will be shifted. The invert flag is set to 0 for fixed shifts, 1 (or anything else) for invertible shifts. Invert flags associated with free shifts will be overwritten inside OCCAM2DMT to point to the locations in the PM vector where the static shifts reside.
Model prejudice file:
The optional prejudice file allows an entire image of the parameter vector to be read, representing a preferred model which the inverted model will try to stay close to. This is accomplished by adding a weighted difference between the inversion model and the preferred model to the penalty measure. Normally the weights for the preferred model are set to zero, but if one wants to regularize into fixed structure, or has other preconceptions about what the resistivity structure should look like, the difference between the model prejudice and the inversion model can be penalized.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only `OCCAM2MTPREJ\_1.0' is supported. Beware of leading blanks after the 18 character label.
Field 2: 18 character `DESCRIPTION' label followed by 80 character description. This can be any information the user wishes to attach to the file. Ignored by OCCAM2DMT.
Field 3: 18 character `NO. PARAMETERS' label followed by number of parameters, NPREJ. This number must match NP in the STARTUP file.
Field 4: NPREJ entries of preferred model parameters, free format.
Field 5: NPREJ entries of integer relative weights to give the preferred model, free format. The weights are multipliers of the default weighting, which is the same as the horizontal weight given to a block. Thus 1 will provide equal weight to the difference between the model and the preferred one as to the roughness (difference between adjacent bricks). 0 will produce no preference for the model prejudice.
Some Notes on Usage
If you haven't done so, read the notes associated with the inversion routines, OCCAMv2.x. Anything said there is applicable to the 2D MT problem.
Since Phil Wannamaker's forward code is at the heart of occam2Dmt, read his notes also. Strictly speaking, his recommendations for the FE mesh should be adhered to. However, smooth inversion puts you at a great advantage over parameterized inversion, in that sharp conductivity boundaries are avoided wherever possible. Also, you can't predict where the conductivity contrasts are going to occur anyway. The net result is that the regularization bricks can be very sparsely gridded. We have done quite well using only two finite elephants per regularization brick width, two in height for the upper two layers followed by one in height below these. The bottom and side bricks have to be meshed a lot better than this, usually with at least half a dozen elephants.
Other Programs
There are a few programs that accompany the academic release of the Occam program and 2DMT subroutines. These are MakeModel2Dmt, which makes a set of generic startup files for a 2D inversion after asking a few questions, FindResponse, which finds the response of a model in an ITERxx file and outputs the model as a set of postscript commands to draw coloured blocks, and SynMod, which generates a set of synthetic data with noise from a startup or ITERxx file. These programs are not very clean, and have no documentation.
Introduction. These routines are based on the original 1D code which has been shoehorned into the new Occam formalism. Because of this they lack a certain neatness, but are fully functional, as far as I know.
Distribution and Portability
The 1D MT inversion needs files `mt1D.f' and `mt1D.inc' in addition to the occam files. Also available are makemodel1dmt.f, which is a stand-alone program to create all the input files you need (except, of course, the data), and plusd.f, Bob Parker's program to find the minimum misfit. One of these days I will modify plusd to take a standard occam data file, but not this month.
The routines are written in FORTRAN 77, and I explicitly take advantage of the FORTRAN 77 feature that states that thou shall not necessarily go through a do-loop.
The common blocks used to convey information from INPUTM() to FORMOD() and OBJMAT() are not declared in the calling program, so you need to make sure that your compiler uses static commons or can be set to do so. Local variables are assumed to be dynamic, and anything that needs to be saved is explicitly SAVEd.
An overview of using occam1Dmt.
a) Set the maximum dimensions for your problem, if they are not already adequate. Dimensions only appear in 'occamdim.inc'.
b) Compile and link occam.f and mt1D.f, making sure that all the include files are available.
c) Create the input files. This will normally be done by another computer program, such as WSE's Geotools or the simple program that I distribute with the academic version of the code. occam1Dmt does only minimal checking to make sure the various inputs are valid and consistent, relying on the code that generates the files to get it right.
There are 2 files required to run the occam1Dmt, in addition to the startup file required by OCCAMv2.0. These are the data file and model file.
Setting the dimensions
The Occam distribution notes describe setting NPP, NDD and NNP in occamdim.inc. Basically, these numbers are the maximum number of model parameters (layers), maximum number of data (frequencies $\times$ components), and maximum number of penalty terms. In the distribution, NNP is set to 2*NPP, which will accommodate vanilla 2D inversion.
The Data File:
This file is pointed to by the startup file, but points to no other files. In the description below, a `field' is usually one line, but not always so.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only `OCCAM1MTDATA\_2.0' is supported. Make sure there are no leading blanks after the 18 character label.
Field 2: 18 character `DESCRIPTION' label followed by 80 character title. This can be any information the user wishes to attach to the file. Ignored by occam1Dmt.
Field 3: 18 character `NUM DATA' label followed by the number of data (free format), NDATA.
Field 4: ` Period Type Datum Std Error' banner. Skipped by OCCAM1DMT.
Field 5: NDATA lines of data (free format). Each line is made up of the period, data type (see below), datum, and standard error. Period is in seconds. The data types are as follows:
The Model File:
This is pointed to by the startup file.
Field 1: 18 character `FORMAT' label followed by the format name. Currently only `OCCAM1MTMOD\_2.0' is supported. Beware of leading blanks aftet the 18 character label.
Field 2: 18 character `MODEL NAME' label followed by 80 character title. This can be any information the user wishes to attach to the file. Ignored by occam1Dmt
Field 3: 18 character `DESCRIPTION' label followed by 80 character model description. This can be any information the user wishes to attach to the file. Ignored by occam1Dmt
Field 4: 18 character `NUM LAYERS' label followed by (free format) number of layers in the model, NLAYER. This number includes (in 1D terminology) the terminating half-space. Layers are counted from the surface down.
Field 5: ` Thickness Resistivity Penalty Prejudice Penalty' banner. Skipped by OCCAM1DMT.
Field 6: NLAYER entries of layer thickness (meters), resistivity, penalty, prejudice, and a second penalty. If resistivity is set less than zero, it is a free parameter for inversion (normal), otherwise it is fixed (I have not yet tested a mixture of fixed an preferred parameters). The first penalty is the weight of the difference between this layer and the one below it (which is why the last entry doesn't matter). The prejudice is a preferred resistivity for this layer, which is applied using the second penalty. The second penalty will normally be zero, so whatever appears in the prejudice column will be ignored.
Document Revision 1.00.
The seafloor1d directory contains all the programs needed to perform a 1D inversion of seafloor data. It works, but is preliminary.
Model file: This is similar in format to the mt1D case, except that the first layer is the water layer, (not counted in the layer count (forth entry in the file)).
Data file: Enties are 4 parameters, datum, and error. Any electric field magnitudes are $V/Am^2$, phases are degrees.
The parameters are: type (see below), freq (Hz), range (m), azimuth (degrees, 0=polar 90=equatorial).
Type refers to the component being inverted: