# Laboratory exercise #2

Document Sample

```					                             Tuesday, April 14, 2009                                                  1

EXPLORATION GEOPHYSICS - MINERAL METHODS
ESS 136B – SPRING 2009
Professor R.L. McPherron

MAGNETIC METHODS - OUTLINE OF LAB #3
Automatic Inversion of a General Dipole
Wednesday, April 15, 2009

1.0 INTRODUCTION
Many years ago there were still vacant areas on the UCLA campus where it was
possible to do geophysics. One of these was at the corner of Gayley and Veteran (SW
corner of UCLA), a site now occupied by the regional dense shelving facility of the UC
library system. Our class performed a gridded survey at this site and discovered a
significant magnetic anomaly. The teaching assistant returned to obtain additional data
on a denser grid. This laboratory exercise is designed to illustrate principles of magnetic
data display and interpretation using MATLAB. Both the data manipulation program and
principles of the magnetic method are important.

2.0 DATA ENTRY AND DISPLAY
The data from the survey were recorded manually on data sheets by the students
performing the survey. This survey was performed on a magnetic N/S - W/E grid at 2
meter spacing. Additional data were acquired later by the teaching assistant along some
lines using a 1-meter grid. Both sets of numbers were entered into a table with black
pencil. To carry out an interpretation additional numbers were interpolated with a hand
calculator using continuity of both the N-S and E-W profiles as a guide. The interpolated
numbers were added to the table in red and shaded. Measured extremes were circled.
Subsequently I interpolated the data to a one meter grid, nngrid2p, with 25 rows (N to S),
and 19 columns (W to E). This two dimensional array has been moved to the class web

http://www.igpp.ucla.edu/public/rmcpherr/Teaching/ESS136B_Spring2009/Labs/Lab03_VacantLotInversion/

2.1 Transfer the Data Array to your Directory

C:\My Teaching\ESS135_136\ESS136B_2009\Labs\Lab03_VacantLotInversion\LAB03_AutoInverseGendip.DOC
October 1, 2012
2

2.2 Create a Publication Quality Contour Map of Data Array
You can significantly improve the quality of the graphs you save with the following code:

fout = ‘Drive:\Path\file.ext’
print(2, '-dmeta', '-r864', fout)

In all of the following keep in mind that the anomaly is at the east edge of the map (and
array) with the negative perturbation to the north (top) and the positive perturbation to the
south (bottom). When you make contour and surface plots make them look like maps
with the correct orientation of the x and y axis. You will have to use the flipud and
fliplr functions to accomplish this.

This is the structure of my program, dolab2.m. Note I have removed some of the code
that I think you should now be able to do leaving the outline and some of the hard parts.
You may use this as a guide. One point of this exercise is to make you think about the
relation between directions on maps and profiles and their relation to the way you store
the data in Matlab.

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%flip the array top to bottom
nn    = flipud(nngrid2p);     %This flips the array top to bottom

   We did the profiles from south to north. We recorded the data points in a map format
with north at the top of the table, west at the left edge. This means that row 1 of the file
is the North edge of the plot. When MATLAB plots a column it starts with row 1 (North)
at left edge of plot. This is backwards to a principal profile (south to north). Also
MATLAB does contour and surface maps with row 1 (North) at the bottom, and row N
at the top of the map. Thus a straight contour of the file nngrid2p would reverse north
and south. Thus we must flip the input top to bottom.

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% MAKE A CONTOUR PLOT
figure(1)

% Use the contourf function
cs = contourf(….);
clabel(cs,'manual')

   If you don’t pass a set of contour levels you won’t get a good contour map. I suggest
you make an array called levs that goes from -200 to 1400 by 50 nT. This plot will
be improved if you also give it the actual values along X- and Y-axis. On the previous
page I tell you the number and spacing of the rows and columns so generate
appropriate vectors for the X (north) and Y (east) locations. (Start at zero, not 1). You
are going to be confused with Matlab indexing and conventions! In magnetic

Page 2
lab95_03.doc
Tuesday, April 14, 2009                                                   3

surveys X is north (up), Y is east (right), and Z is down. In Matlab X is horizontal
(right) and Y is vertical (up), and Z axis (up) is used to display contours of height of
the function of x and y. What this means is that you must reverse the positions of x
and y in the calls to these functions.
   Label some of the contour labels with the clabel function. Use help so that you
understand these two commands.

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%MAKE A SURFACE PLOT
figure(2)
Use the surfc function

   The surfc function plots the data properly with the origin of x and y coordinates at
the bottom center of the image with east to the right and up and north to the left and
up. The problem is that you can’t see the negative part of the anomaly because it is
hidden. You can rotate the graph with the view command. I used a 180 rotation
viewed from 35 degrees above placing the origin at the back of the plot. Use “help
view”

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%MAKE PLOT OF THE PRINCIPLE PROFILE
figure(3)
Use
Use the regular plot command and extract the column with largest
signature

   Examine the contour map (or the original array) and determine which column is
directly over the source of the anomaly. This is the principal profile. Put this column
in an array pp = ???? and them plot it as a function of x (north distance).
%ADD A REGIONAL TREND USING GINPUT AND POLYFIT
disp('Click left and right edges for linear trend')
[l,r] = ginput(2);        %Points at left and right on linear trend
cof   = polyfit(l,r,1);   %Coefficients of linear fit
% Evaluate this trend at every measurement point using polyval
fit   = polyval(….);      %A fit to all measurement points
Add this fit to the plot

   I made the south end of profile a little above the trend and the north end on the trend.
%SUBTRACT THE FITTED LINEAR TREND TO OBTAIN MEASURED ANOMALY
anom = ????

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%MAKE A PLOT OF THE MEASURED ANOMALY
figure(4)

%DETERMINE THE EXTREMA LOCATIONS AND VALUES
[ppx,ix]= max(anom);
[ppn,in]= min(anom);
xex   = xe(ix);       %North coord of maximum

C:\My Teaching\ESS135_136\ESS136B_2009\Labs\Lab03_VacantLotInversion\LAB03_AutoInverseGendip.DOC
October 1, 2012
4

xen       = xe(in);       %North location of minimum
hw        = abs(xex-xen);     %Separation of extrema (~depth)
x0        = mean([xex,xen]); %Approximate location of dipole

S{1} = sprintf('xmax=%4.1f, Bmax=%5.0f',xex,ppx);
S{2} = sprintf('xmin=%4.1f, Bmon=%5.0f',xen,ppn);
S{3} = sprintf('x0=%4.1f, hw=%4.1f',x0,hw);
text(15,600,S,'FontSize',11,'FontWeight','bold')

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%ITERATE THE FORWARD MODEL TO GET AN APPROXIMATE FIT
   This anomaly seems to have a magnetic north south strike from the rough symmetry
around the N-S direction. But it can not be an induced dipole because the negative
part of the anomaly is too great. Assume that this is a buried dipole under the
principal profile and use the function gendip.m to model it.
   The assumption of symmetry means muy is zero.
   The fact that the positive peak is to the south suggests the dipole is pointed down like
an induced dipole would be implying muz is positive.
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%ITERATE THE FORWARD MODEL TO GET AN APPROXIMATE FIT
% Create the spatial array defining the principal profile
xp = [xe',zeros(size(xe'))];        %Create 2-col array of surf loc

% Set the fixed parameters
global fixparm vr
fixparm = ???;      Fixed parameters (inc=60, comp = 0 {total field})
vr      = ???;                  %Normalization parms set to all 1.0

% Make initial guess of the variable parameters
mux = ???;
muy = 0;
muz = ???;                              %Start values for dipole
b   = [x0,0,hw,mux,muy,muz];            %First guess of variable
parameters, Note I used the half width as depth guess.

Loop to do the forward model iteration
lop = 1;
hold off
while(lop==1)
dbi      = gendip(b,xp);       %Model with first guess parms
pef      = 1- var(anom-dbi)/var(anom);
plot(xp( ,1),anom,'b-o',xp( ,1),dbi,'r-x')
fprintf(1,' %5.1f%% %5.1f %5.1f %5.1f ... %5.0f %5.0f %5.0f\n',100*pef,b)
fprintf(1,'Enter a new “b” array and CR to proceed or lop=0 and CR to
end\n')
keyboard
End

You terminate this loop by setting lop = 0.

Use the MATLAB help command to study the command "contour". Contour displays an
array with row 1 at the BOTTOM of screen and column 1 at the left. Use the flipud

Page 4
lab95_03.doc
Tuesday, April 14, 2009                                                5

command to reverse the array top to bottom so it contours properly as a map. Construct
three arrays to input as arguments to the contour command. First create distance along
the vertical screen axis (the y axis to MATLAB). In our geophysical coordinate system
this is distance x towards magnetic north. The array xn = [0:1:24]; has 25 elements. Next
create distance along the horizontal screen axis (x axis in MATLAB) the geophysical y or
east axis. The array ye = [0:1:18]; has 19 elements. Finally, create the array of contour
levels between -200 and +1400 nT. Use a 100 nT contour level, levs = [-200:100:1400];.
Now contour the data using the arrays just created. Remember to reverse xn, and ye in
the command argument string as x and y geophysical are rotated 90 degrees relative to the
screen. Add informative labels, title, annotation, and print out your final plot.

2.2 Make A Principle South To North Profile Plot

Extract the single column “jj” from nngrid2p containing the maximum and
minimum magnetic anomalies. Also select the row “kk” with the peak.
pp = nngri2p(:,jj);       %This is the N-S principal profile
bpp = flipud(pp);         %Make it south to north
bew = nngridp(kk,:);      %This is the west to east profile
Plot these profiles on a graph with appropriate labels and titles. Add a graph of the west
to east profile through the maximum disturbance. Annotate the two traces on the graph.

2.Do a Regional Trend Analysis on Principal Profile

Use “ginput” to select a good linear trend in the S-N direction. Add the fit to the
plot of the principal profile. Annotate traces, label axes, add title, plot, print your plot. I
actually did this by using ginput, polyfit, and polyval on the principal profile. See my
code above.

2.4 Plot the Principle Anomaly Profile

I decided that the bets fit regional trend for the entire map was a constant 580 nT.
Subtract this number from every point on your map to get the anomaly.

3.0 FUNCTION GENDIP

I have provided this function to you in today’s handout at the website. Remember the
“global” command and get the parameters in the correct order.

3.3 Coding a Main Program to do Forward Models
Write a simple loop to try different values for mux and muz. Each time through show the
data and the model. Calculate the prediction efficiency. I have provided most of this code
above.

C:\My Teaching\ESS135_136\ESS136B_2009\Labs\Lab03_VacantLotInversion\LAB03_AutoInverseGendip.DOC
October 1, 2012
6

4.0 AUTOMATIC INVERSION OF PRINCIPLE PROFILE
Use nlinfit with your final forward model to get a better solution. Make a plot of this as
well.
%Nonlinear inversion of the principal profile
beta = nlinfit(xp,anom,'gendip',b);

Use the output of your manual forward model as the initial guess to this.

Note that Rachel discovered that the nlinfit in your version of Matlab (2008b) does not
converge. I don’t know why! My version (2008a) does converge. I have placed my
version in the class directory. You must download this function and place it in your
directory where your script for this lab is located.

5.0 EVALUATE MODEL FOR OPTIMUM PARAMETERS AND COMPARE
DATA AND MODEL IN A PLOT
Plot the anomaly determined above in blue with circle symbols. Add the plot of your optimum model
obtained just above. Calculate the prediction efficiency of you model.

WRITEUP

Title Page - Modify your previous one

Abstract - Concise summary of entire report

Introduction - Rephrase what I have said above as Introduction and add any thing you
remember of what I said in class. Be brief.

Body of Report - Same outline as this handout. Here I expect the description to be an
extended Figure caption for each figure I requested. These can be terse. Just tell me
what is in the figures.

Conclusions - This should be connected text and state in plain English what you set out
to do, what you did, what you found, and what it means. This is the conclusion
section for your first report to your employer who sent you out to find a toxic waste
dump.

Appendices - Attachments that I gave you and listings of your programs.

Page 6
lab95_03.doc
Tuesday, April 14, 2009                                    7

6. APPENDICES TO HANDOUT FOR LAB EXERCISE #3, SPRING 1995
nngrid2p.asc
24    588   521    452   387    352   370    416      462    498    523      540   554   564     569   573        582   590   592   590
23    639   559    473   386    336   352    405      453    490    516      535   550   561     567   572        581   590   596   600
22    695   605    509   411    352   367    417      458    487    511      530   544   554     560   566        576   588   598   607
21    752   660    566   473    418   427    462      481    491    506      522   532   537     541   548        563   580   595   609
20    777   694    613   539    491   482    493      499    499    496      494   495   498     503   515        536   560   582   602
19    743   680    622   571    525   485    469      488    500    471      430   415   420     431   452        486   523   555   584
18    676   639    605   574    530   468    432      455    473    423      342   284   260     269   312        382   456   512   556
17    613   597    579   557    522   469    427      417    406    353      253   137     21     8      83       220   362   459   524
16    574   555    535   514    491   466    434      391    348    309      222    67    -98   -140    -60       116   302   423   500
15    567   513    469   444    436   438    427      388    350    334      294   203     80    21      53       182   334   434   497
14    564   483    427   413    422   429    425      410    401    414      440   461   456     415   384        406   454   492   520
13    542   479    450   474    501   484    454      452    476    522      610   747   860     877   817        718   628   582   566
12    525   508    523   581    616   574    517      517    562    630      751   942   1118   1187   1127       955   771   666   614
11    539   568    611   670    697   656    603      603    647    716      825   978   1119   1185   1142       988   813   706   646
10    566   615    665   711    729   705    675      676    711    772      849   932   1000   1031   1003       906   793   714   661
9    583   618    655   693    715   711    701      709    741    794      852   896   929     919   896        832   765   710   664
8    588   609    633   661    683   692    697      713    745    794      843   874   879     858   822        782   742   701   660
7    587   623    649   656    658   668    684      705    736    782      827   851   849     821   783        752   723   691   655
6    583   632    659   653    641   649    670      692    722    765      806   828   825     797   761        733   708   680   650
5    578   603    620   624    627   639    657      678    706    744      782   806   807     782   747        719   695   670   644
4    573   569    574   593    616   632    645      663    688    720      753   779   786     765   732        705   682   660   638
3    566   562    568   586    608   621    631      646    666    690      718   744   754     738   709        685   665   648   631
2    560   570    580   591    602   611    620      631    644    661      681   701   709     696   673        656   643   632   622
1    556   570    582   590    597   606    615      621    627    636      647   654   653     638   621        615   614   613   611
0    553   564    574   583    592   604    613      615    614    615      614   606   590     571   560        567   582   592   599
0     1     2      3     4      5     6      7        8      9      10       11    12     13    14      15        16    17    18   19

C:\My Teaching\ESS135_136\ESS136B_2009\Labs\Lab03_VacantLotInversion\LAB03_AutoInverseGendip.DOC
October 1, 2012
8

gendip.m
function dbi = gendip(b,xp)
% Modified April 15, 2003 to use with nlinfit.m
%    b       =       varparms     = [x0,y0,d,mux,muy,muz]
%    This program has a global statement and expects the calling program
to
%    have the same statement
% function dbi=gendip(xp,fp,vp,vr)
%_______________________________________________________________________
%|        dbi = gendip(xp,fp,vp,vr)
%| Feb 15, 1993 1900 PST      Mods: April 20, 1995 RLM
%|m-file to calculate magnetic effect of general dipole
%| xp =      normalized col array (N x 2] of obs vectors [x,y]
%| fp =      fixed parameters [incb, comp] inclination of the B field and
a
%            switch comp = 1 =>'Z' field and comp = 0 => 'B'
%| vp =      normalized variable parameters [x0,y0,d,mux,muy,muz]
%| vr =      parameters used to normalize the vp parameters
%| dbi =     Anomaly in chosen component
%|......Normalized distance is input, ie xp/xr
%|       Normalized variable parameters are input, ie x0/xr, ..., mux/mr
....
%|......NOTE! The partial derivitive routine expects only one output
variable
%| You must change this routine for Z or F output
%|______________________________________________________________________
_

global fp vr                    %Duplicate this in the calling program

comp   = fp(2);                 %Get the switch
dr = pi/180;

% DISASSEMBLE THE FIXED PARAMETERS
incb    = fp(1);

% DISASSEMBLE THE NORMALIZATION FACTORS
xr = vr(1);
mr = vr(4);                         %B inclin, dist. scale, Mu scale

% DISASSEMBLE THE VARIABLE PARAMETERS
x0 = b(1);
y0 = b(2);
d   = b(3);                     %First three elements of variable vec
mux = b(4);
muy = b(5);
muz = b(6);                 %Second three elements of variable vec
%
xd = (xp(:,1)-x0)/d;
yd = (xp(:,2)-y0)/d;
l   = sqrt(xd.^2+yd.^2+1);
con = (mr/xr^3)*(100/d^3);      %First term compensates for scaling
if(comp==0)
%......Code for total field
ci = cos(dr*incb); si=sin(dr*incb);
ter5    = 3*(mux*xd+muy*yd-muz).*(ci*xd-si)./l.^5;

Page 8
lab95_03.doc
Tuesday, April 14, 2009                                                9

ter3    = (mux*ci+muz*si)./l.^3;
both    = ter5-ter3;
dbi = con*both;
elseif(comp==1)
%......Code for Z component
ser5    = 3*(mux*xd+muy*yd-muz)./l.^5;
ser3    = (muz)./l.^3;
zoth    = -(ser5+ser3);
dbi= con*zoth;
else
disp('Unknown comp! Expect 0 fo B or 1 for Z - STOP')
end
%disp('gendip'),keyboard
return

C:\My Teaching\ESS135_136\ESS136B_2009\Labs\Lab03_VacantLotInversion\LAB03_AutoInverseGendip.DOC

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 7 posted: 10/2/2012 language: Unknown pages: 9
How are you planning on using Docstoc?