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

球体磁异常正演计算程序matlab磁异常正演是地球物理勘探中的一种重要方法，它通过计算地下物质对地球磁场的扰动，来推断地下物质的性质和分布。

球体磁异常正演是磁异常正演中的一种常见方法，它假设地下物质为球体形状，通过计算球体对地球磁场的扰动，来推断球体的磁性和位置。

在球体磁异常正演计算中，matlab是一种常用的计算工具。

下面介绍一种基于matlab的球体磁异常正演计算程序。

我们需要定义球体的参数，包括球心坐标、半径和磁化强度。

假设球心坐标为(x0,y0,z0)，半径为R，磁化强度为M，我们可以用matlab的向量表示：x0 = 0;y0 = 0;z0 = 1000;R = 500;M = 0.1;接下来，我们需要定义计算区域的参数，包括计算网格的大小和范围。

假设计算区域为一个正方形，边长为L，我们可以用matlab的向量表示：L = 5000;nx = 101;ny = 101;x = linspace(-L/2,L/2,nx);y = linspace(-L/2,L/2,ny);[X,Y] = meshgrid(x,y);然后，我们需要定义地球磁场的参数，包括地球磁场的强度和方向。

假设地球磁场的强度为B0，方向为沿z轴正方向，我们可以用matlab的向量表示：B0 = 50000;Bx = zeros(nx,ny);By = zeros(nx,ny);Bz = B0*ones(nx,ny);接下来，我们可以用球体磁异常正演的公式计算球体对地球磁场的扰动。

假设我们采用了文献[1]中的公式，我们可以用matlab的代码表示：r = sqrt((X-x0).^2 + (Y-y0).^2 + (z0).^2);cos_theta = z0./r;sin_theta = sqrt(1-cos_theta.^2);cos_phi = (X-x0)./sqrt((X-x0).^2 + (Y-y0).^2);sin_phi = (Y-y0)./sqrt((X-x0).^2 + (Y-y0).^2);Br = (M/4/pi)*((3*cos_theta.^2-1)./r.^3);Btheta = (M/4/pi)*(3*sin_theta.*cos_theta./r.^3);Bphi = (M/4/pi)*(-3*sin_phi.*cos_theta./r.^3);我们可以将球体对地球磁场的扰动与地球磁场相加，得到总磁场。

TECHNOLOGY AND INFORMATION80 科学与信息化2023年10月下基于Python和Fortran程序实现瞬变电磁法一维反演软件的设计与开发孙大利 李方舟 费靖淇中国地震局第一监测中心 天津 300171摘 要 瞬变电磁方法可以有效地识别地下水分布情况，从而为地下城市空间科学安全的开发利用提供有力的帮助。

为了降低瞬变电磁法数据处理分析的使用门槛，使非专业人员能够快速分析瞬变电磁探测数据，本文基于Python和Fortran语言实现瞬变电磁法中心回线系统一维Occam反演软件核心算法和交互界面的设计开发，简化了专业处理流程，更高效地发挥了瞬变电磁法探测在城市地下空间开发中的辅助作用。

关键词 Occam反演；应用软件；地下空间Design and Development of One-Dimensional Inversion Software for Implementing Transient Electromagnetic Method Based on Python and Fortran Programs Sun Da-li, Li Fang-zhou, Fei Jing-qiThe First Monitoring and Application Center of China Earthquake Administration, Tianjin 300171, ChinaAbstract Transient electromagnetic method can effectively identify the distribution of groundwater, so as to provide powerful help for the scientific and safe development and utilization of underground urban space. In order to reduce the threshold of use of transient electromagnetic data processing and analysis, and enable non-professionals to quickly analyze transient electromagnetic detection data, this paper implements the design and development of the core algorithm and interactive interface of the one-dimensional Occam inversion software for implementing transient electromagnetic method central loop system based on Python and Fortran language, simplifies the professional processing process, and plays the auxiliary role of transient electromagnetic detection in the development of urban underground space more efficiently.Key words Occam inversion; application software; underground space引言近年来随着城市化的不断发展，城市规模的不断扩大，地上有限的空间资源难以满足人们日常生活、工作的需求。

Matlab电磁时间反演（EMTR）进行空间-变压器局部放电定位EMTR法原理可以参照以下这个程序来进行书写% 变电站PD定位仿真（TR法）(4个传感器)close all;clear all;clc;f0=500e6; %中心频率v=3e8; %声速t_low=0; t_up=200e-9; N=1001; %时间范围if mod(N,2)==0;N=N+1;endt=linspace(t_low,t_up,N)';Lx=20; %变压器长Ly=20; %变压器宽Lz=20; %变压器高derta_x=0.3; %传感器间隔,排列方式：从左至右，从上至下x_coor=-0.5*derta_x+derta_x*(0:1);z_coor=zeros(1,2);y_coor=-0.5*derta_x+derta_x*(0:1);x=repmat(x_coor,[1,2]);z=repmat(z_coor,[1,2]);y=[y_coor(1)*ones(1,2),y_coor(2)*ones(1,2)];x0=-7; y0=5; z0=5; %PD位置t_lag1=150e-9; %参考传感器延迟时间（传感器1）w=20e-9; %高斯脉冲脉宽A=1; %高斯脉冲幅度t_lag=t_lag1+[sqrt((x-x0).^2+(y-y0).^2+(z-z0).^2)-sqrt((x(1)-x0).^2+(y(1)-y0).^2+(z(1)-z0).^2)]/v; %各个传感器相对于传感器1的延时x01=10;y01=-5;z01=3;t_lag3=120e-9;t_lag2=t_lag1+[sqrt((x-x01).^2+(y-y01).^2+(z-z01).^2)-sqrt((x(1)-x01).^2+(y(1)-y01).^2+(z(1)-z01).^2)]/v; %各个传感器相对于传感器1的延时% x01=10;y01=-5;z01=3;pd=A*exp(-4*log(2)*(t-t_lag).^2/(w^2)).*sin(2*pi*f0*(t-t_lag))+A*exp(-4*log(2)*(t-t_lag2).^2/(w^2)).*sin(2*pi*f0*(t- t_lag2)); %各个传感器PD信号% pd=A*exp(-4*log(2)*(t-t_lag).^2/(w^2)); %各个传感器PD信号(宽带信号效果不好)% pd1=A*exp(-4*log(2)*(t-t_lag1).^2/(w^2)); %传感器1PD信号% pd2=A*exp(-4*log(2)*(t-t_lag2).^2/(w^2)); %传感器2PD信号% pd3=A*exp(-4*log(2)*(t-t_lag3).^2/(w^2)); %传感器3PD信号% pd4=A*exp(-4*log(2)*(t-t_lag4).^2/(w^2)); %传感器4PD信号% pd4=0;pd_reversed=fliplr(pd')'; %时间反演N_distance1=200;N_distance2=200;theta=0:pi/N_distance1:pi/2;phy=0:pi/N_distance2:pi*2;M=[];f=(0:(N-1)/2)'/(diff(t(1:2)))/N;n=0; progress=0;for j=1:length(theta)for k=1:length(phy)t_delay=(x*cos(theta(j))*cos(phy(k))+y*cos(theta(j))*sin(phy( k))+z*sin(t heta(j)))/v;t_delay=t_delay-t_delay(1);T=exp(1i*2*pi*f*t_delay);T((N+1)/2+1:N,:)=conj(T((N+1)/2:-1:2,:));pd_test=real(ifft(fft(pd_reversed).*T));M(j,k)=sum((sum(pd_test')).^2);% M=[M,sum(pd_test.^2)];n=n+1;endendif N_distance1~=1M=(M-min(min(M)))/(max(max(M))-min(min(M)));subplot(1,2,1)surf(phy/pi*180,theta/pi*180,M);xlabel('方位角/度')ylabel('仰角/度')zlabel('归一化')shading interp;colormap('jet');% subplot(2,1,2)%contour(linspace(0,360,N_distance),linspace(0,90,N_distance),M) end% 局放源定点[~,YI]=max(max(M));phy_max=phy(YI);[~,YI]=max(max(M'));theta_max=theta(YI);N_distance3=500;if (phy_max>pi/2)&&(phy_max<3*pi/2)dx=linspace(0,-Lx/2,N_distance3);elsedx=linspace(0,Lx/2,N_distance3);enddy=dx*tan(phy_max);dz=sqrt(dx.^2+dy.^2)*tan(theta_max);n1=length(find(dx<=Lx/2));n2=length(find(dy<=Ly/2));n3=length(find(dz<=Lz));N_dis=min([n1,n2,n3]);dx=dx(1:N_dis);dy=dy(1:N_dis);dz=dz(1:N_dis);M1=[];for j=1:N_dist_delay=sqrt((dx(j)-x).^2+(dy(j)-y).^2+(dz(j)-z).^2)/v; t_delay=t_delay-t_delay(1);T=exp(-1i*2*pi*f*t_delay);T((N+1)/2+1:N,:)=conj(T((N+1)/2:-1:2,:));pd_test=real(ifft(fft(pd_reversed).*T));M1(j)=sum((sum(pd_test')).^2);end[~,YI]=max(M1);subplot(1,2,2)plot3(dx,dy,dz)hold onplot3(dx(YI),dy(YI),dz(YI),'bp','Markersize',10) xlabel('x/米') ylabel('y/米')zlabel('z/米')plot3(x0,y0,z0,'r*')plot3(x0,y0,z0,'ro')。

基于MATLAB的地震正演模型实现贾跃玮;杨锐【摘要】人工合成地震正演模型是进行三维模型计算的基础,在地震勘探领域具有重要意义.针对地震勘探原理,作者运用MATLAB强大数学计算和图像可视化功能,对一个三层介质模型制作二维人工合成地震记录.文章说明地震记录形成的物理机制,介绍了地质模型的构造及参数选择,针对具体地质模型制作合成地震记录,成功验证了褶积模型原理.【期刊名称】《计算机与数字工程》【年(卷),期】2009(037)007【总页数】4页(P132-135)【关键词】地震;MATLAB;正演【作者】贾跃玮;杨锐【作者单位】中国地质大学地下信息探测技术与仪器教育部重点实验室,北京,100083;川庆钻探工程有限公司地质勘探开发研究院,成都,610059【正文语种】中文【中图分类】工业技术【文献来源】https:///academic-journal-cn_computer-digital-engineering_thesis/0201241737639.html总第 237 期2009 年第 7 期计算机与数字工程 Computer&DigitalEngineering Vol.37No.7132基于 1\/IATLAB的地震正演模型实现+贾跃玮¨ 杨 ~ 孙（中国地质大学地下信息探测技术与仪器教育部重点实验室” 北京100083 ）（川庆钻探工程有限公司地质勘探开发研究院” 成都 610059 ）摘要人工合成地震正演模型是进行三维模型计算的基础，在地震勘探领域具有重要意义。

针对地震勘探原理，作者运用 MATLAB强大数学计算和图像可视化功能，对一个三层介质模型制作二维人工合成地震记录。

文章说明地震记录形成的物理机制，介绍了地质模型的构造及参数选择，针对具体地质模型制作合成地震记录，成功验证了褶积模型原理。

关键词地震 MATI 。

AB正演中图分类号P631.4ASeismicForward ModelingBasedonMATLAB Jia Yuewei" YangRuiz) (Geo-detectionLab,Ministry of EducationChinaUniversityof Geosciences".Beijing 100083) (GeologicalExplorationandDevelopmentResearchInstituteof SichuanPetroleumAdministration2) ,Chengdu 610059) Abstract A syntheticseismicforwardmodelis thebaseof3Dmodelingcalculation,andhasprofoundmeaningin seismicexploration.Basedontheprinciplesofseismicexploration,theauthorappliedgreatmathematiccalculationandfigural visualization of MATLABtobuild asyntheticseismic recordingof athree-layermodel.Thephysicalprincipleofseismic explorationisintroducedatthe firstpart of the article.Then,thetectonic structureof the geophysicalmodelandparametersaredelineated.In thelast partof thepaper,the methodof syntheticseismic recordsbasedonthe specifiedge-ophysicalmodelispresented. Keywordsseismic,MATLAB,forward modelingClassNumber P631.4／／ 1 引言r地震勘探是利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应，推断地下岩层的性质和形态的地球物理方法。

大地电磁测深数据一维自动反演作者：赵长海来源：《价值工程》2011年第25期摘要：在学科实践中，反演是一种相当普遍的手段。

地球物理反演则是通过对地面、地下、空间，甚至海洋上的观测（地震仪、重力仪、地电及地磁仪）资料进行分析计算，来推断地球内部介质的地震波速度、密度、电导率等参数的分布，从而得到地球内部介质分布的二维或三维结构图像。

本文根据阮百尧教授在其“电阻率激发极化法测深数据的一维最优化反演方法”一文中所述，导出了在大地电磁测深中与电阻率测深一维最优反演相似的反演方法，通过自动迭代反演出地下模型参数。

计算表明，这是一种快速和实用的自动反演方法。

Abstract: In practice, inversion is a method used widely. Geophysics inversion deduces the distribution of seismic-wave speed, density, conductivity and so on in earth interior, then obtains the 2D or 3D structure images for the distribution of earth interior medium, by analyzing the observed date on surface, in subsurface and in space, even on the sea. Referring to the correlate articles by Profession Ruan B.Y, Derived in the magnetotelluric sounding and resistivity sounding similar to the one-dimensional inversion method optimal inversion, Inverted through the automatic iterative model parameters. Calculations show that this is a fast and practical method for the automatic inversion.关键词：大地电磁测深；水平层状；一维正反演；快速反演Key words: magnetotelluric sounding；horizontal layered；1D forward and inversion；rapid inversion中图分类号：P631.3文献标识码：A 文章编号：1006-4311（2011）25-0285-020 引言对于简单的地电模型（一维）一般能够计算其解析解。

大地电磁一维正反演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('******正演参数不够******！')endu(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_endT1=linspace(10^t,10^(t+1),Num_DT);for i=1:1:Num_DTj=j+1;T(j)=T1(i);endendfor j=1:1:(T_end-T_start+1)*Num_DTfor 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);endrou_T(j)=rou(1)*(u(1)^2+v(1)^2);endfigure(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));endfclose(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);function dif(x)[T,rou]=mt1d(5,[5,10,15,10,5],[10,20,30,100],-3,6,20);log_rou=log10(rou);log_T=log10(T);log_omeg=log10(2*3.14159265./T);for i=2:1:length(T)-1d_log_rou(i)=log_rou(i+1)-log_rou(i-1);d_log_T(i)=log_T(i+1)-log_T(i-1);d_log_omeg(i)=log_omeg(i+1)-log_omeg(i-1);a(i)=(d_log_rou(i))/(d_log_T(i));fai(i)=(d_log_rou(i))/(d_log_omeg(i));enda(1)=(log_rou(2)-log_rou(1))/(log_T(2)-log_T(1));a(length(T))=(log_rou(length(T))-log_rou(length(T)-1))/(log_T(length(T))-log_T(length(T)-1)); fai(1)=(log_rou(2)-log_rou(1))/(log_omeg(2)-log_omeg(1));fai(length(T))=(log_rou(length(T))-log_rou(length(T)-1))/(log_omeg(length(T))-log_omeg(length( T)-1));bsdk_rou=rou.*((1+a)./(1-a));bsdk_d=rou.*T;bsdk_d=bsdk_d/(8*(pi^2)*10^(-7));bsdk_d=sqrt(bsdk_d);bsdk_d=sqrt((rou.*T)./(8*(pi^2)*10^(-7)));s=bsdk_d./rou;hold onplot(log10(T),bsdk_rou,'--ks','MarkerFaceColor','k','MarkerSize',10);figure(2)plot(log10(bsdk_rou),log10(bsdk_d),'--ks','MarkerFaceColor','k','MarkerSize',10);。

实验报告课程名称：地电学课题名称：大地电磁层状模型数值模拟实验专业：地球物理学姓名：xx班级：06xxxx完成日期：2016 年11月26日目录一、实验名称 (3)二、实验目的 (3)三、实验要求 (3)四、实验原理 (3)五、实验题目 (4)六、实验步骤 (4)七、实验整体流程图 (8)八、程序及运行结果 (9)九、实验结果分析及体会 (14)一、实验名称大地电磁层状模型数值模拟实验二、实验目的（1）学习使用Matlab编程，并设计大地电磁层状模型一层，二层，三层正演程序（2）在设计正演程序的基础上实现编程模拟（3）MATLAB软件基本操作和演示.三、实验要求（1）利用MT一维测深法及其相关公式，计算地面上的pc视电阻率和ph相位，绘制视电阻率正演曲线和相位曲线并分析。

（2）利用Matlab软件作为来实现该实验。

四、实验原理（一）、正演的概念：正演是反演的前提。

在实际地球物理勘探中，一些模型的参数是不容易确定的，如埋藏在地下的地质体模型的岩性、厚度、产状等参数，我们把这些描述未知模型的参数的集合定义为“模型空间”。

为了获得这些模型参数，可以利用那些可以直接观测的量来推测，而这些能够直接观测的量的集合则被称作“数据空间”。

如果把模型空间中的一个点定义为m，把数据空间中的一个点定义为d，按照物理定律，可以把两者的关系写成式中，G为模型空间到数据空间的一个映射。

我们把给定模型m求解数据d的过程称为正演问题。

（二）、MT一维正演模型简介大地电磁法作为一种电磁类勘探方法，它的模型参数为一组能够表征地球物理勘探目标体的电性参数，即目标体电阻率和相应层的层厚度。

所谓一维模型，即介质在三维空间中沿两个方向上模型参数是不变的，只在另一个方向上特征属性会变化。

在此一维模型即指水平层状一维介质，即介质只在沿垂直于地面上的方向上电性（电阻率）变化，在另外两个方向上保持不变的典型特征，所以就构成一组电阻率不同的电性层，抽象出来即是一组由电阻率及对应的层厚度构成的电性层数。

大地电磁测深一维正反演摘 要 本文推导了大地电磁测深的理论计算表达式，并以水平层状介质为例，利用推导的正演计算式在MATLAB 软件平台上进行正演，比较了不同层介质参数的视电阻率曲线。

简要介绍了阻尼最小二乘法反演的基本原理和反演迭代步骤，并对多种层介质进行了反演。

关键词 大地电磁，一维正反演，阻尼最小二乘法1 引 言20世纪50年代初，苏联学者吉洪诺夫和法国学者卡尼亚的经典著作奠定了大地电磁测深法（MT ）的基础。

它是利用大地仲频率范围很宽（4410~10-Hz ）广泛分布的天然变化的电磁场，进行深部地质构造研究的一种频率域电磁测深法。

由于该法不需要人工建立场源，装备轻便、成本低，且具有比人工源频率测深法更大的勘探深度，所以除主要用于研究地壳和上地幔地质构造外，也常被用来进行油气勘查、地热勘探以及地震预报等研究工作。

几十年来，由于大地电磁测深法具有以下几个优点：不受高阻屏蔽，对低阻分辨率高；不用人工供电，勘探成本低且工作方便；勘探深度范围大。

使大地电磁法在矿产勘探及普查、地壳岩石圈电性结构研究、海洋地球物理勘探、地热勘探、能源勘探、隐伏岩溶水结构、天然地震预测等都扮演着至关重要的角色。

大地电磁也存在一些缺点，比如在实际应用的过程中整理后的数据存在分散的情况；频率范围不够宽，特别是缺少高频成分，受噪音影响大信噪比低；所需观察时间长，致使野外工作效率低。

随着基础理论、技术手段、仪器设备的不断完善和发展，进一步改进和解决这些问题，才能将大地电磁法更好的应用于生产服务当中。

2 视电阻率及水平地层大地电磁测深曲线的理论计算方法 2.1大地电磁测深理论的几点假设和论证吉洪诺夫和卡尼亚提出了假设并论证了以下几点：①将场源近似地看为平面电磁波垂直入射大地。

②引入波阻抗的概念（Z=E/H ），表征地球电性分布对大地电磁场的响应。

③利用单点大地电磁场观测研究地球电性分布是可能的。

2.2视电阻率及水平地层上的理论计算表达式视电阻率概念是从均匀介质中电阻率和波阻抗关系引申出来的。