一、均匀半空间模型 1) h1→∞，可得ρa=ρ1，cth(ik1h1)→1 说明均匀半空间所得视电阻率就是真电阻率，是 一个常数，即ρ1，因此后边结果中将不再对均匀半空间情况做讨论。 2) 而相位也是常数（-45°），这也说明了模型是电性均一的，相位差不变。 二、.两层地电模型 1) 视电阻率公式为： clc; clear; close all; n=2; p=[100,500]; h=[100,0]; f=logspace(-4,4,20); F=2*pi*f; u=4*pi*10^(-7); L=length(f); R=linspace(0,0,n); ps=linspace(0,0,L); k=linspace(0,0,n-1); i=sqrt(-1); xw=linspace(0,0,L); for j=1:L; R(n)=1; %层数 %P 为层电阻率 %磁导率 %事先定义相位 for a=n-1:-1:1 %第 a 层复波数 k(a)=sqrt(-2i*pi*f(j)*u/p(a)); t1=sqrt(p(a+1))*R(a+1)+sqrt(p(a))*(1-exp(-2i*k(a)*h(a)))/(1+exp(-2i*k(a)*h(a))); t2=sqrt(p(a))+sqrt(p(a+1))*R(a+1)*(1-exp(-2i*k(a)*h(a)))/(1+exp(-2i*k(a)*h(a))); R(a)=t1/t2; ps(j)=p(1)*(abs(R(1)))^2; end ps(j)=p(1)*(abs(R(1)))^2; ps(j)=ps(j)/p(1); F(j)=sqrt(10^(7)*p(1)/f(j))/h(1); Z1=-sqrt(i*2*pi*f(j)*u*p(1)); xw(j)=atan(imag(Z1*R(1))/real(Z1*R(1)))*180/pi; end %x=min(F):0.001:max(F); %y=spline(F,wx,x);%做插值，光滑曲线 loglog(F,xw,'*',F,xw,'b'); title('大地电磁电测深波阻抗相位曲线'); xlabel('频率(HZ)'); ylabel('波阻抗相位'); grid on; 

%暂停 box on; pause %x=min(F):0.001:max(F); %y=spline(F,ps,x);%做插值，光滑曲线 semilogx(F, ps,'*',F, ps,'c'); title('大地电磁电测深视电阻率曲线'); xlabel('频率(HZ)'); ylabel('视电阻率'); grid on; box on; %磁导率 function [ps,xw]=MT_1(rho,h) u=(4e-7)*pi; T=logspace(-5,3,30); F=2*pi./T; i=sqrt(-1); k=zeros(size(rho,2),size(T,2)); for N=1:size(rho,2) k(N,:)=sqrt(-2*i*pi*u./(T.*rho(N))); %k(N,:)=sqrt(-2*i*pi*u./(T.*rho(N))-u*8.8419e-012*omega.^2) end m=size(rho,2); y=-(i*u*2*pi)./(T.*k(m,:)); for nn=m-1:-1:1; A=-(i*2*pi*u)./(T.*k(nn,:));B=exp(-2*k(nn,:)*h(nn)); y=A.*(A.*(1-B)+y.*(1+B))./(A.*(1+B)+y.*(1-B)); end ps=(T./(u*2*pi)).*(abs(y).^2); xw=-atan(imag(y)./real(y)).*180/pi; b=[ps,xw]; %x=min(T):0.001:max(T); %y=spline(T,wx,x);%做插值，光滑曲线 loglog(T,xw,'- rx'); title('大地电磁电测深波阻抗相位曲线'); xlabel('频率(HZ)'); ylabel('波阻抗相位'); grid on; 

%暂停 box on; pause semilogx(T, ps,'- rx'); title('视电阻率曲线图'); xlabel('频率(HZ)'); ylabel('视电阻率Ω.m'); grid on; box on;