实验2--磁性体磁场正演程序
磁场实验的测量步骤与技巧
磁场实验的测量步骤与技巧磁场实验是物理学中常见的实验之一，通过测量磁场的强度和方向，可以帮助我们理解磁场的特性和应用。
在进行磁场实验时，我们需要遵循一定的测量步骤并掌握一些技巧，以确保实验结果的准确性和可靠性。
本文将介绍磁场实验的测量步骤和一些实用的技巧。
一、准备实验器材和环境在进行磁场实验之前，首先需要准备好实验所需的器材和环境。
常用的实验器材有磁铁、磁针、霍尔效应传感器等。
确保这些器材的品质良好，并定期检查和校准。
同时，为了减小外界磁场对实验的干扰，应在无明显磁场干扰的环境中进行实验，避免有大型电流或磁性材料的附近。
二、测量磁场强度测量磁场强度是磁场实验的重要内容。
以下是一般的测量步骤：1. 将磁场测量仪器和传感器放置在待测点附近，确保传感器的位置准确。
可以利用三维坐标系来标定传感器的位置，从而确定测量点的坐标。
2. 打开磁场测量仪器，校准传感器，使其与磁场方向垂直。
校准时，可以利用已知磁场强度的标准磁铁或磁场源进行调整。
3. 进行实际测量。
可以通过改变传感器与被测磁场的相对位置或者调节磁场强度来进行实验。
记录测量值，并考虑实验误差。
4. 重复实验多次，取平均值，以提高结果的准确性。
三、测量磁场方向除了测量磁场强度，我们也可以测量磁场方向。
以下是一般的测量步骤：1. 准备一个磁针。
磁针应具有一定的灵敏度，可以指示磁场的方向。
2. 将磁针静置一段时间，使其与周围环境达到热平衡。
确保磁针完全静止。
3. 将磁针放置在待测点附近。
注意磁针的位置应该与磁场方向垂直。
4. 观察磁针的指向，记录磁场的方向。
磁针会指向磁场的南北极。
5. 若需要测量平面内的磁场方向，可以旋转磁针的位置，并观察磁针的指向变化。
通过多次观察，可以确定磁场的方向。
四、注意事项及技巧在进行磁场实验时，还需要注意以下事项，并掌握一些实用的技巧：1. 避免外界磁场干扰。
在实验中，应尽可能减小外界磁场的干扰。
例如，远离大型电流和磁性材料。
磁性物理实验指导书
磁性物理实验讲义磁性物理课程组编写电子科技大学微电子与固体电子学院二O一二年九月目录一、起始磁导率温度特性测量和居里温度测试计算分析 (1)二、电阻率测试及磁损耗响应特性分析 (3)三、磁致伸缩系数测量与分析 (6)四、磁化强度测量与分析 (9)五、磁滞回线和饱和磁感应强度测量 (11)六、磁畴结构分析表征 (12)一、起始磁导率温度特性测量和居里温度测试计算分析(一) 、实验目的：了解磁性材料的起始磁导率的测量原理，学会测量材料的起始磁导率，并能够从自发磁化起源机制来分析温度和离子占位对材料起始磁导率和磁化强度的影响。
(二）、实验原理及方法：一个被磁化的环型试样，当径向宽度比较大时，磁通将集中在内半径附近的区域分布较密，而在外半径附近处，磁通密度较小，因此，实际磁路的有效截面积要小于环型试样的实际截面。
为了使环型试样的磁路计算更符合实际情况，引入有效尺寸参数。
有效尺寸参数为：有效平均半径r e ，有效磁路长度l e ，有效横截面积A e ，有效体积V e 。
矩形截面的环型试样及其有效尺寸参数计算公式如下。
⎪⎪⎭⎫ ⎝⎛-=211211ln r r r r r e （1）⎪⎪⎭⎫ ⎝⎛-=211211ln2r r r r l e π （2）⎪⎪⎭⎫ ⎝⎛-=2112211ln r r r r h A e （3）e e e l A V = （4）其中：r 1为环型磁芯的内半径，r 2为环型磁芯的外半径，h 为磁芯高度。
利用磁芯的有效尺寸可以提高测量的精确性，尤其是试样尺寸不能满足均匀磁化条件时，应用等效尺寸参数计算磁性参数更合乎实际结果。
材料的起始磁导率（i μ）可通过对环型磁心施加线圈后测量其电感量（L ）而计算得到。
计算公式如式（5）所示。
20i e eA N L l μμ=（5）其中：μ0为真空磁导率，4π×10-7 H·m -1；N 为线圈匝数。
磁性材料起始磁导率（µi ）的定义式如式（6）所示。
二、磁性体磁场正演
M s = ( M x2 + M z2 )1/ 2 = M (cos 2 I cos 2 A′ + sin 2 I ) −1 M z −1 is = tg = tg (tgI sec A′) Mx
四、程序编写与图示结果
例子：1 例子：1）先定义测区范围
% 测点分布范围 dx=5; % X方向测点间距 dy=5; % Y方向测点间距 nx=81; % X方向测点数 ny=81; % Y方向测点数 xmin=-200; % X方向起点 ymin=-200; % Y方向起点 x=xmin:dx:(xmin+(nx-1)*dx); % X方向范围 y=ymin:dy:(ymin+(ny-1)*dy); % Y方向范围 [X,Y]=meshgrid(x,y); % 转化为排列
四、程序编写与图示结果
例子：4）计算磁异常（Za、 例子：4）计算磁异常（Za、….)
hax=(u*m*((2*x.^2-y.^2-h.^2)*cos(I)*cos(a)3*h*x.*sin(I)+3*x.*y.*cos(I)*sin(a)))./(4*pi*(x.^2+y.^2+h.^2).^(5/2)); hay=(u*m*((2*y.^2-x.^2-h.^2)*cos(I)*sin(a)3*h*y.*sin(I)+3*x.*y.*cos(I)*cos(a)))./(4*pi*(x.^2+y.^2+h.^2).^(5/2)); za=(u*m*((2*h.^2-x.^2-y.^2)*sin(I)-3*h*x.*cos(I)*cos(a)3*h*y*cos(I)*sin(a)))./(4*pi*(x.^2+y.^2+h.^2).^(5/2)); dt=hax.*cos(I)*cos(a)+hay.*cos(I)*sin(a)+za.*sin(I);%总磁场强度异常
磁性体磁场正演
§3、规则形体的磁场
薄板状体
薄板状体可看作是厚板的特殊 情况。在磁法中“厚”与“薄”也 是一个相对概念。在一定限度 内当板状体的b<<h 时,称其 为薄板，反之为厚板。 厚板与薄板的剖面曲线形态类 似。薄板的磁场表达式可从厚 板的磁场表达式简化导出。 厚板状体可以看作薄板状体组 合而成，薄板的异常窄，幅值 小，而厚板异常宽，幅值大。
H ax
μ 0 M s • sin α rB ln = 2π rA
μ 0 M s • sin α Za = (Δϕ ) 2π
§3、规则形体的磁场
倾斜磁化板状体磁场
斜磁化指板的侧面与磁化强 度Ms斜交的情况，γ≠0 斜交磁化厚板的顶面、底面 和侧面都要出现磁荷。 斜交磁化无限延伸厚板磁场 Za图形随磁化倾角：
x = 0, Z a⊥ = Z a max H a⊥ = 0 μ0 2ms = 4π R 2
规则形体的磁场
四、水平圆柱体
通常将自然界中延深和宽度都比较小，沿走向很长 的磁性体看作水平圆柱体。 一.水平圆柱体的磁场表达式: 若为垂直磁化,即is=90°,或I=90°时：
μ0 ms ( R 2 − x 2 ) Z a⊥ = 2π ( x 2 + R 2 ) 2
磁性体磁场正演
规则形体的磁场
球体的重力异常：Δg = GM
h (x + h )
2 2 3 2
规则形体的磁场
球体的重力异常：Δg = GM
h (x + h )
2 2 3 2
磁异常垂直分量 Z a
qm h ： Z a = 4π 2 2 3 2 (x + h )
规则形体的磁场
一、单极的磁场（顺轴磁化、无限延深柱体）
5第五章 磁性体的磁场
(5.2-4)
• 对二度体(即走向为无限长的物体)，因 引力位与坐标变量y无关，故有
0 H ax M xVxx M zVxz 4 G
Hay 0
(5.2-5)
0 Za M xVxz M zVzz 4 G
2 0 2 a
1/2
• 对上式两端取平方，并除以 B
2 2
2 0
，则得
(5.2-15)
T T Ba Ba 2 2 cos B0 B0 B0 B0
• 当Bａ＜＜Ｂ0时，上式中的平方项可略去。 • 例如，在中纬度地区，B０＝50000ｎＴ，若 Bａ≤2000ｎＴ时，则(Bａ/Ｂ０）２ ≤0.0016。 • 又因Δ Ｔ＜Ｂa，故(Δ Ｔ／Ｂ０）２项也可 以略去。 • 因此，(5.2-15)可以简化为
uu r uu r uu r M 1 1 M gradQ divQ divQ M r r r
• 将上式代入(5.2-8)式，得
dv M 1 U divQ dv divQ M 4 v r v r
第二节 计算磁性体磁场的基本公式
• 计算磁性体磁场的方法有多种，主要是体积分公式法、 面积分公式法和重磁位的泊松公式法。 • 一般应根据形体特征，选用简便的基本公式。 • 在讨论各种磁性体的磁场时，我们选用右手直角坐标， XOY面为观测平面； • Y轴沿地质体走向方向，X轴垂直走向，为测线方向，Z 轴垂直向下，坐标原点选在磁性体中心或顶面中心在 地 面的投影点。 • 我们还规定磁性体以外的计算点以P表示，该点坐标以 x、y、z表示； • 磁体内的体元点以Q表示，其坐标为ξ 、ε 、δ 。
磁异常的处理、解释及应用
上面，简单概述了磁异常各类正问题及其正演方法。其中，均匀磁化规则形体正（演） 问题、正演方法及场的解析表达式，是磁法的基础，具有重要的理论意义和实际意义将重点 讨论。
第二节 有效磁化强度矢量与总磁场异常 Δ T 的一般表达
一、有效磁化强度矢量
已知总磁化强度矢量由感应磁化强度与剩余磁化强度两矢量组成。设总磁化强度矢量 M
K 这种磁性体的参数 k 和 M r 需用张量来描述，其正演问题是磁法中最复杂的正问题。从
70 年代后期，国内外学者相继研究出一些数值正演方法。我国学者把有限元和边界元等数 值计算方法引用到这一复杂正演问题中来，取得了一系列有理论和实际价值的成果。
（五）磁场的模拟测定
前述各类正问题的求解还可以通过实验室模拟测定来解决。模拟测定方法分为静磁场 模拟方法与低频交变场模拟方法。实践已经证明，两类模拟测定方法是可行的。
3、组合体近似法：把磁化强度均匀或分区均匀的任意形态磁性体，用多个均匀磁化规 则形体的组合形体近似代替；各个均匀磁化规则形体的磁化强度可以相同或不同。该磁性体 磁场的近似值，等于各规则体解析场值之和。作为组件的规则形体有正方体、直立长方体、 倾斜长方体、有限长水平n棱柱体等。因为直立长方体的多个ln项可以合并成一项计算，而 且在一定条件下多个tg-1项亦可合并计算，使计算速度大大加快，又因其组合任意形体的能 力较强，故直立长方体组合法得到了普遍应用。
基于MATLAB的磁性体磁场正演模拟研究
1.引言 根 据 已 知 地 质 体 及 磁 性 体 的 形 态 、质 量 及 磁 性 、空间等
分 布 来 计 算 其 磁 场 分 布 的 过 程 ，称 为 磁 场 正 演 问 题 。磁场正
在 进 行 磁 性 体 正 演 模 拟 时 ，通常会设置一定的假设条 件 ，如磁性体为简单规则形体、磁性体被均勻磁化、只研究
下转第104页
上接第102页
⑵真空过減机组增加喷淋 真空过滤机组过滤硫磺流程示意图如图3所 示 ，吸收氧 化塔塔底硫磺浆通过硫磺浆泵输送到真空过滤机，滤液通过 滤液泵返回至吸收氧化塔。
厂；研究方向：螯合铁天然气脱硫工艺和克劳斯法天然气脱硫工艺 装置投产运行研究。
104 C当he代nm化ica工lIn研ter究mediate
技术应用与研究
2 0 1 8 •0 8
单 个 磁 性 体 、观 测 面 水 平 、磁 性 体 的 剩 余 磁 化 强 度 为 零 。本
-7
x 10
文 以 球 体 为 例 进 行 磁 场 正 演 模 拟 ，具 体 情 况 如 下 ：
参数 、走 向 等 ），掌握磁性体磁场正演的计算方法，
： 中图分类号 T
： 文献标识码 A
Research on Magnetic Field Forwarding Modeling of Magnetic Body Based MATLAB
演 是 磁 场 反 演 的 基 础 ，野 外 测 量 的 最 后 成 果 是 磁 异 常 的 等 直 线 平 面 图 和 剖 面 图 ，可 以 反 映 地 下 目 标 物 与 围 岩 的 磁 性 差 异 所 产 生 的 磁 异 常 特 征 。磁 勘 探 就 是 利 用 磁 异 常 反 演 求 地 下 目 标物的磁性参量和几何参量及其空间分布位置。
磁法勘探 05 磁异常的正演
∆T = Ta cosθ = Ta cos(Ta , T0 ) &
4
上式表明，当磁异常 强度不大时 可近似把△ 看作是 强度不大时， 看作是Ta在 上式表明，当磁异常Ta强度不大时，可近似把△T看作是 在T0方向的 投影；航空磁测中一般Ta<1000nT，在进行高精度地面磁测的地区，一般 投影；航空磁测中一般 ，在进行高精度地面磁测的地区，一般Ta 也不大。因此将△ 近似看作 近似看作Ta在 方向的投影， 足够的精度。另外， 也不大。因此将△T近似看作 在T0方向的投影，有足够的精度。另外，T0 相当大的区域内，方向是不变的（ 内变化1º左右），因此 左右），因此， 在相当大的区域内，方向是不变的（10000km2内变化 左右），因此，可以 看作是Ta在固定方向的投影 的物理意义与Za、 类似 类似， 把△T看作是 在固定方向的投影。这样△T的物理意义与 、Ha类似，都 看作是 在固定方向的投影。这样△ 的物理意义与 在固定方向的分量。 是Ta在固定方向的分量。 在固定方向的分量
1 V = Gσ ∫∫∫ dv r v
U= −1 M ⋅ grad p 4π
∫∫∫
v
1 dv r
G为引力常数， 为密度，将引力位公式带入磁位公式，可得： 为引力常数， 为密度，将引力位公式带入磁位公式 可得： 磁位公式， 为引力常数 σ
U =− 1 M ⋅ grad p V 4πGσ
µ0 [M xVxx + M yVyx + M zVzx ] 4πGσ µ0 [M xVxy + M yVyy + M zVzy ] H ay = 4πGσ µ0 [M xVxz + M yVyz + M zVzz ] Za = 4πGσ
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《应用地磁学》实验报告姓名：张嘉琪学号： 1010112225指导教师：李淑玲实验地点：实验室319实验日期： 2014-05-24实验二：磁性体磁场正演一、实验目的：1、通过球体、水平圆柱体磁场的正演计算，掌握简单规则磁性体正演磁场的计算方法；2、通过计算认识球体与水平圆柱体磁场的一般分布规律，了解影响磁性体磁场的主要因素（如磁性体的形体、物性参数、走向或计算剖面的选择等），培养学生实际动手能力与分析问题的能力。
二、实验内容用Matlab 语言或C 语言编程实现球体和水平圆柱体的磁场(包括Za 、Ha 、Δt)的正演计算。
三、实验要求假设地磁场方向与磁性体磁化强度方向一致且均匀磁化的情况下,当地磁场T=50000nT ，磁倾角I=60°，球体与水平圆柱体中心埋深R=30m ，半径r=10m ，磁化率k=0.2（SI ），计算（观测）剖面磁化强度水平投影夹角A ′=0°时：1、正演计算球体的磁场（Za 、Hax 、Hay 、ΔT ），画出对应的平面等值线图、曲面图及主剖面异常图；2、正演计算水平圆柱体的磁场（Za 、Ha 、ΔT ），画出主剖面异常结果图；3、通过改变球体与水平圆柱体的几何参数、磁化强度方向（I ）、计算剖面的方位角(A ′)，观察主剖面磁场Za 的变化，分析磁化方向与计算剖面对磁性体磁场特征的影响。
四、实验原理球体与水平圆柱体磁场（Za 、Ha 、ΔT ）的计算公式是以磁化强度倾角I 、有效磁化倾角is 和剖面与磁化强度水平投影夹角A ′来表达。
1、球体磁场的正演公式： [[[⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎪⎬⎫'-'---++='+-'--++='+-'--⋅++=]sin cos 3cos cos 3sin )2( )(4]cos cos 3sin 3sin cos )2( )(4]sin cos 3sin 3cos cos )2( )(42222/522202222/522202222/52220A I Ry A I Rx I y x R R y x m Z A I xy I Ry A I R x y R y x m H A I xy I Rx A I R y x R y x m H a ay ax πμπμπμ()]sin 2sin 32sin cos 3cos 2sin 3sin cos )2(cos cos )2(sin )2[(42222222222222222/52220A I yR A I xy A I xR A I R x y A I R y x I y x R R y x mT '-'+'-'--+'--+--++=∆πμ2、水平圆柱体磁场的正演公式：⎪⎪⎭⎪⎪⎬⎫+-+-=--+=]sin 2cos )[()(12]cos 2sin )[()(12222220222220s s s a s s s a i Rx i x R R x m H i Rx i x R R x m Z πμπμ()()()()[]902cos 2902sin sin sin 2222220----+=∆s s s si Rx i x R i I R x m T πμ 3、有效磁化强度Ms 与有效磁化倾角is ：⎪⎭⎪⎬⎫'==+'=+=--)sec ()sin cos (cos )(112222/122A tgI tg M M tg i I A I M M M M x z s z x s 五、计算程序代码：1、球体matlab 代码：clc;clear;%% 测点分布范围dx=5; % X 方向测点间距dy=5; % Y 方向测点间距nx=81; % X 方向测点数ny=81; % Y 方向测点数xmin=-200; % X 方向起点ymin=-200; % Y 方向起点x=xmin:dx:(xmin+(nx-1)*dx); % X 方向范围y=ymin:dy:(ymin+(ny-1)*dy); % Y 方向范围[X,Y]=meshgrid(x,y); % 转化为排列% 球体参数i=pi/3; %磁化倾角ia=0; %剖面磁方位角R=10; % 球体半径 mv=4/3*pi*R^3u=4*pi*10^(-7);%磁导率T=0.5*10^(-4);%地磁场强度k=0.2;%磁化率M=k*T/u; %磁化强度 A/mm=M*v; %磁矩D=30; % 球体埋深 m% 球体Za理论磁异常Za=(u*m*((2*D.^2-X.^2-Y.^2)*sin(i)-3*D*X.*cos(i)*cos(a)-3*D*Y.*cos(i)*sin(a))). /(4*pi*(X.^2+Y.^2+D.^2).^(5/2));% 球体Hax理论磁异常Hax=(u*m*((2*X.^2-Y.^2-D.^2)*cos(i)*cos(a)-3*D*X.*sin(i)+3*X.*Y.*cos(i)*sin(a)) )./(4*pi*(X.^2+Y.^2+D.^2).^(5/2));%球体Hay理论磁异常Hay=(u*m*((2*Y.^2-X.^2-D.^2)*cos(i)*sin(a)-3*D*Y.*sin(i)+3*X.*Y.*cos(i)*cos(a)) )./(4*pi*(X.^2+Y.^2+D.^2).^(5/2));%球体ΔT理论异常T=Hax*cos(i)*cos(a)+Hay*cos(i)*sin(a)+Za*sin(i);%绘平面异常等值线图(二维)figure(1),clf,subplot(221),contourf(X,Y,Hax);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体Hax异常');axis equal,axis([-50 50 -50 50]),colorbar;subplot(222),contourf(X,Y,Hay);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体Hay异常');axis equal,axis([-50 50 -50 50]),colorbar;subplot(223),contourf(X,Y,Za);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体Za异常');axis equal,axis([-50 50 -50 50]),colorbar;subplot(224),contourf(X,Y,T);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体ΔT异常');axis equal,axis([-50 50 -50 50]),colorbar;%绘制曲面图（三维）figure(2),clf,subplot(221),mesh(X,Y,Hax),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Hax异常'),colorbar; subplot(222),mesh(X,Y,Hay),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Hay异常'),colorbar; subplot(223),mesh(X,Y,Za),shading interp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Za异常'),colorbar;subplot(224),surf(X,Y,T),shading interp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体ΔT异常'),colorbar;%绘制主剖面异常等值线Za1=(u*m*((2*D.^2-x.^2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.^2+D.^2).^(5/2) );Hax1=(u*m*((2*x.^2-D.^2)*cos(i)*cos(a)-3*D*x.*sin(i)))./(4*pi*(x.^2+D.^2).^(5/2 ));Hay1=(u*m*((-x.^2-D.^2)*cos(i)*sin(a)))./(4*pi*(x.^2+D.^2).^(5/2));T1=Hax1*cos(i)*cos(a)+Hay1*cos(i)*sin(a)+Za1*sin(i);figure(3),clf,subplot(221)plot(x,Za1,'g-','linewidth',1.3);xlabel('X(m)'),ylabel('理论球体Za异常'); subplot(222)plot(x,Hax1,'k-','linewidth',1.3);xlabel('X(m)'),ylabel('理论球体Hax异常'); subplot(223)plot(x,Hay1,'r-','linewidth',1.3);xlabel('X(m)'),ylabel('理论球体Hay异常'); subplot(224)plot(x,T1,'b-','linewidth',1.3);xlabel('X(m)'),ylabel('理论球体ΔT异常');%绘制异常剖面图figure(4),clf,for i=0:pi/6:pi/2Za2=(u*m*((2*D.^2-x.^2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.^2+D.^2).^(5/2) );hold onplot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(磁倾角改变)'),grid on;endh=legend('Za');legend(h,'boxoff');figure(5),clf,for a=0:pi/6:piA=pi/3;Za2=(u*m*((2*D.^2-x.^2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.^2+D.^2).^(5/2) );hold onplot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(磁方位改变)'),grid on;endh=legend('Za');legend(h,'boxoff');figure(6),clf,for i=pi/3;a=0;R=10:5:20v=4/3*pi*R^3m=M*v;Za2=(u*m*((2*D.^2-x.^2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.^2+D.^2).^(5/2) );hold onplot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(球体半径)'),grid on;endh=legend('Za');legend(h,'boxoff');圆柱体程序代码：clc;clear;%% 测点分布范围dx=5; % X方向测点间距dy=5; % Y方向测点间距nx=81; % X方向测点数ny=81; % Y方向测点数xmin=-200; % X方向起点ymin=-200; % Y方向起点x=xmin:dx:(xmin+(nx-1)*dx); % X方向范围y=ymin:dy:(ymin+(ny-1)*dy); % Y方向范围[X,Y]=meshgrid(x,y); % 转化为排列% 水平圆柱体参数i=pi/3; %磁化倾角a=0;%剖面磁方位角Is=(tan(tan(i)*sec(a)))^(-1);R=10; % 圆柱体横截面半径 mS=pi*R^2; %圆柱体横截面面积u=4*pi*10^(-7); %磁导率T=0.5*10^(-4);%地磁场强度k=0.2;%磁化率M=k*T/u; %磁化强度 A/mMs=M*((cos(i)*cos(a))^2+(sin(i))^2);m=Ms*S; %单位长度的有效磁矩D=30; % 圆柱体中心点埋深 m% 圆柱体Za理论磁异常Za=(u*m*((D.^2-X.^2)*sin(Is)-2*D*X.*cos(Is)))./(2*pi*(X.^2+D.^2)^2);% 圆柱体Ha理论磁异常Ha=(-u*m*((D.^2-X.^2)*cos(Is)+2*D*X.*sin(Is)))./(2*pi*(X.^2+D.^2)^2);%圆柱体ΔT理论异常T=(u*m*sin(i)*((D.^2-X.^2)*cos(2*i-pi)-2*D*X.*sin(2*Is-pi/2)))./(sin(Is)*((D.^2 -X.^2)*sin(2*Is-pi/2)-2*D*X.*cos(2*Is-pi/2)));%绘平面异常等值线图(二维)figure(1),clf,subplot(221),contourf(X,Y,Za);xlabel('X(m)'),ylabel('Y(m)'),title('理论圆柱体Za异常');axis equal,axis([-200 200 -200 200]),colorbar;subplot(222),contourf(X,Y,Ha);xlabel('X(m)'),ylabel('Y(m)'),title('理论圆柱体Ha异常');axis equal,axis([-200 200 -200 200]),colorbar;subplot(223),contourf(X,Y,T);xlabel('X(m)'),ylabel('Y(m)'),title('理论圆柱体ΔT异常');axis equal,axis([-200 200 -200 200]),colorbar;%绘制曲面图（三维）figure(2),clf,%clf清除图形subplot(221),mesh(X,Y,Za),shading interp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论圆柱体Za异常'),colorbar;subplot(222),mesh(X,Y,Ha),shading interp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论圆柱体Ha异常'),colorbar;subplot(223),surf(X,Y,T),shading interp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论圆柱体ΔT异常'),colorbar;%主剖面视图figure(3),clf,subplot(311)for x=-200:5:200Za1=(u*m*((D.^2-x.^2)*sin(Is)-2*D*x.*cos(Is)))./(2*pi*(x.^2+D.^2)^2);hold on;plot(x,Za1,'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体Za异常'); endsubplot(312)for x=-200:5:200Ha1=(-u*m*((D^2-x^2)*cos(Is)+2*D*x*sin(Is)))/(2*pi*(x^2+D^2)^2);hold on;plot(x,Ha1,'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体Ha异常'); endsubplot(313)for x=-200:5:200T1=(u*m*sin(i)*((D^2-x^2)*cos(2*i-pi)-2*D*x*sin(2*Is-pi/2)))/(sin(Is)*((D^2-x^2 )*sin(2*Is-pi/2)-2*D*x*cos(2*Is-pi/2)));hold on;plot(x,T1,'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体T异常');end%绘制异常剖面图figure(4),clf,for i=0:pi/6:pi/2Is=(tan(tan(i)*sec(a)))^(-1);Ms=M*((cos(i)*cos(a))^2+(sin(i))^2);m=Ms*S;Za1=(u*m*((D.^2-X.^2)*sin(Is)-2*D*X.*cos(Is)))./(2*pi*(X.^2+D.^2)^2); hold onplot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常(磁倾角)'),grid on; endh=legend('Za');legend(h,'boxoff');figure(5),clf,for a=0:pi/6:pi/2i=pi/3;Is=(tan(tan(i)*sec(a)))^(-1);Ms=M*((cos(i)*cos(a))^2+(sin(i))^2);m=Ms*S;Za1=(u*m*((D.^2-X.^2)*sin(Is)-2*D*X.*cos(Is)))./(2*pi*(X.^2+D.^2)^2); hold onplot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常（方位角）'),grid on; endh=legend('Za');legend(h,'boxoff');figure(6),clf,for R=10:5:20i=pi/3;a=0;S=pi*R^2;m=Ms*S;Za1=(u*m*((D.^2-X.^2)*sin(Is)-2*D*X.*cos(Is)))./(2*pi*(X.^2+D.^2)^2); hold onplot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常（半径）'),grid on; endh=legend('Za');legend(h,'boxoff');六、实验结果：球体实验结果：平面等值线图：曲面图：主剖面视图：球体参数改变后的主剖面视图：（半径改变）磁倾角改变后的主剖面视图：剖面方位角改变后的剖面图：圆柱体实验结果：平面等值线图：曲面图：主剖面视图：球体参数改变后的主剖面视图：（半径改变）磁倾角改变后的主剖面视图：剖面方位角改变的异常图：七、结果分析：球体分析：平面图特征：球体的磁场ΔT不仅与地磁场方向有关，还与观测剖面有关。