MT1111D_FWD.rar_大地电磁一维_大地电磁一维正演_正演_电磁正演

%此程序为一维层状各向同性介质正演程序，结果得到视电阻率和相位 function[]=MT1111D_FWD(rho,h) %输入参数rho为各层电阻率值，h为各层厚度值 mu=(4e-7)*pi; %常量 %T=logspace(-3,4,40); %产生T为从10的-3次方到10的4之间按对数等分的40个元素的行向量 %f=[320 265 229 194 159 132 115 97 79 66 57 49 40 33 27.5 22.5 18.8 16.2 13.7 11.2 9.4 8.1 6.9 5.6 4.7 4.1 3.4 2.81 2.34 2.03 1.72 1.41 1.17 1.02 0.86 0.7 0.59 0.51 0.43 0.35] %f=10.^(-1:0.1:3) %ff=log(f) f=logspace(-2,2,80); T=1./f; k=zeros(size(rho,2),size(T,2)); %产生k为行数为电阻率rho个数，列数为T个数（在这里为40） %的全0矩阵 for N=1:size(rho,2) %for循环求出从1至N层（N即为电阻率rho的个数）中每一层关于不同T值的传播常数k k(N,:)=sqrt(-i*2*pi*mu./(T.*rho(N))); %在这个for循环 end m=size(rho,2); %m等于rho的个数即层数 Z=-(i*mu*2*pi)./(T.*k(m,:)); %最下面一层不同T时的特征阻抗 for n=m-1:-1:1 %for循环为地面阻抗递推公式 A=-(i*2*pi*mu)./(T.*k(n,:)); B=exp(-2*k(n,:)*h(n)); Z=A.*(A.*(1-B)+Z.*(1+B))./(A.*(1+B)+Z.*(1-B)); end rho_a=(T./(mu*2*pi)).*(abs(Z).^2) %视电阻率 phase=-atan(imag(Z)./real(Z)).*180./pi %相位 figure(11) semilogx(f,rho_a,'k-') %axis([0 1000 0 200]) legend('理论值'); xlabel('f/Hz'); ylabel('\rho(\Omega\cdotm)'); %imwrite(f,'figure1.tif','resolution',[6000,6000]); figure(12) semilogx(f,phase,'k-') %axis([0 1000 0 90]) legend('理论值'); xlabel('f/Hz'); ylabel('phase(\circ)'); %imwrite(f,'figure2.tif','resolution',[6000,6000]);