【MATLAB】球体磁异常正演

最新推荐文章于 2024-04-18 10:12:04 发布

PourRevenir

最新推荐文章于 2024-04-18 10:12:04 发布

阅读量942

点赞数 3

文章标签：数学建模 matlab

本文链接：https://blog.csdn.net/2201_75598815/article/details/130850457

版权

利用场源体在上半空间的磁场表达式计算其磁场值，这里简单讨论球体的磁异常计算。

球体磁场的正演公式：

$H_{ax}=\frac{\mu_{0}}{4\pi} \frac{m}{(x^{2}+y^{2}+R^{2})^{5/2}}[(2x^{2}-y^{2}-R^{2})cosIcosA'-3RxsinI+3xycosIsinA']$

$H_{ay}=\frac{\mu_{0}}{4\pi} \frac{m}{(x^{2}+y^{2}+R^{2})^{5/2}}[(2y^{2}-x^{2}-R^{2})cosIcosA'-3RysinI+3xycosIsinA']$

$H_{ax}=\frac{\mu_{0}}{4\pi} \frac{m}{(x^{2}+y^{2}+R^{2})^{5/2}}[(2h^{2}-x^{2}-y^{2})cosIcosA'-3RxsinI-3RycosIsinA']$

$\Delta T = H_{ax}cosIcosA'+H_{ay}cosIsinA'+Z_{a}sinI$

式中R为埋深；m为球体磁矩，且m = MV（M为磁化强度，V为球体体积）；I为磁倾角；A‘为磁偏角。

计算程序如下：

clear;

% 声明参数
k = 0.015;                  % 磁化率
T = 50000;                  % 磁场强度（nT）
R = 50;                     % 中心埋深（m）
r = 30;                     % 球体半径（m）

u0 = 4 * pi * (10e-7);      % 真空磁导率
V = (3 * pi * r.^3) / 4;    % 球体体积
M = k * T;                  % 球体磁化强度
m = V * M;                  % 球体磁矩

% 绘图范围
x = -200:1:200;             %测点
y = -200:1:200;             %测线
[X,Y] = meshgrid(x,y);

open = 1;

while open

sc = '输出图像：磁倾角I 磁偏角A';
disp(sc)

I = input('I= ');
A = input('A= ');

% 计算磁异常
Hax = u0.*m.*((2*X.^2-Y.^2-R.^2).*cosd(I).*cosd(A)-3.*R.*X.*cosd(I).*sind(A)+3.*X.*Y.*cosd(I).*sind(A))./(4.*pi.*((x.^2+y.^2+R.^2).^(2/5)));
Hay = u0.*m.*((2*Y.^2-X.^2-R.^2).*cosd(I).*cosd(A)-3.*R.*Y.*cosd(I).*sind(A)+3.*X.*Y.*cosd(I).*cosd(A))./(4.*pi.*((x.^2+y.^2+R.^2).^(2/5)));
Za  = u0.*m.*((2*R.^2-X.^2-Y.^2).*sind(I)-3.*R.*X.*cosd(I).*sind(A)+3.*R.*Y.*cosd(I).*sind(A))./(4.*pi.*((x.^2+y.^2+R.^2).^(2/5))); 

% 绘制图像
figure('Name','球体磁异常');
subplot(2,2,1);
surfc(Hax);shading flat
colormap(jet)
xlabel('X/m')
ylabel('Y/m')
zlabel('Hax')
title('水平面磁异常Hax');

subplot(2,2,3);
surfc(Hay);shading flat
xlabel('X/m');
ylabel('Y/m');
zlabel('Hay');
title('水平面磁异常Hay');

subplot(2,2,4);
surfc(Za);shading flat
xlabel('X/m');
ylabel('Y/m');
zlabel('Za');
title('垂直面磁异常Za');

open = input('再次运行输入1,结束输入0: ');
end

I=40，A’=60时计算结果如下：