大地电磁一维正反演MATLAB程序

大地电磁一维正反演MATLAB程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

% Purpose:MagnetoTelluric one dimensional forward modeling %

% Ming-Cai ZHang 17st,OCt,2008 CSU-IPGE % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

% variable declaration %

% Input: %

% 1、n---->>>The number of layer %

% 2、rou-->>>Density for every layer %

% 3、h---->>>Thickness for every layer %

% 4、T_start---->>>start time %

% 5、T_end----->>>end time %

% 6、Num_DT-->>the number of sampling time in time interval %

% Output: %

% 1、rou_T-->>apparent resistivity at time % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

function [T,rou_T]=mt1d(n,rou,h,T_start,T_end,Num_DT) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

% Control subjacent variable to build model %

%n=2; %

%rou(1:n)=[200,600]; %

%h(1:(n-1))=10; %

% Control subjacent variable to change continued time %

%T_start=-3; %

%T_end=4; %

%Num_DT=5; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

np=nargin;

if (np<=5)

error('******正演参数不够******！')

end

u(1:n)=1;

v(1:n)=0;

j=0;

T=zeros(1,(T_end-T_start+1)*Num_DT);

for t=T_start:1:T_end

T1=linspace(10^t,10^(t+1),Num_DT);

for i=1:1:Num_DT

j=j+1;

T(j)=T1(i);

end

end

for j=1:1:(T_end-T_start+1)*Num_DT

for m=(n-1):-1:1

[P(m+1),Q(m+1)]=set_pqm1(rou(m+1),rou(m),u(m+1),v(m+1));

G(m)=(4*3.14159265*h(m))/sqrt(10*rou(m)*T(j));

e_G(m)=exp(-G(m));

[A(m),B(m)]=set_ABm(P(m+1),Q(m+1),G(m),e_G(m));

u(m)=(1-(A(m))^2-(B(m))^2)/((1+A(m))^2+B(m)^2);

v(m)=-2*B(m)/((1+A(m))^2+B(m)^2);

end

rou_T(j)=rou(1)*(u(1)^2+v(1)^2);

end

figure(1)

picture1=plot(log10(T),(rou_T),'--*r','MarkerFaceColor','k','MarkerSize',10); %picture1=plot(log10(T),(rou_T),'--rs','LineWidth',2,...

%'MarkerEdgeColor','k',...

%'MarkerFaceColor','b',...

%'MarkerSize',5);

warndlg('请保存正演数据，方便以后利用！！');

pathname='e:\mt1d';

outfile1=fullfile(pathname,'forward.txt');

fid=fopen(outfile1,'w');

fprintf(fid, '%s\n','T Log(T) Resitivity');

for i=1:1:length(T)

fprintf(fid,'%6.2f %12.8f %12.8f\n',T(i),log10(T(i)),rou_T(i));

end

fclose(fid);

function [Am,Bm]=set_ABm(Pm1,Qm1,Gm,e_Gm)

Am=e_Gm*(Pm1*cos(Gm)-Qm1*sin(Gm));

Bm=e_Gm*(Qm1*cos(Gm)+Pm1*sin(Gm));

function [pm1,qm1]=set_pqm1(roum1,roum,um1,vm1)

a1=(sqrt(roum1/roum)*um1);

a=a1^2;

b1=(sqrt(roum1/roum)*vm1);

b=b1^2;

pm1=(1-a-b)/((1+a1)^2+b);

qm1=(-2*b1)/((1+a1)^2+b);