圆锥体完全均衡下重力异常正演 [MATLAB]

  在完全均衡的模型下,若地表有一圆锥体(山峰等),计算跨越山顶的截面上所得到的各种重力异常。

 

地壳密度 $kg\cdot m^{-3}$上地幔密度 $g\cdot cm^{-3}$地表地形圆锥体半径 (km)地表地形圆锥体高度 (km)计算莫霍面形变圆锥半径 (km)计算莫霍面形变圆锥高度 (km)地壳厚度 (km)
$2.8\times 10^{3}$$3.5\times 10^{3}$$2.0$$2.0$$2.0$$8.0$$30.0$

 

  计算结果如下。横坐标单位:m,纵坐标单位:mGal

 

重力异常

  MATLAB代码如下:

% 生成地下体的布格重力异常
syms r a z x h;
f = r / sqrt((z + h)^2 + r^2 + x^2 - 2*r*x*cos(a))^3;

dense = - 700; G = 6.67E-11;
depth = 30000; subheight = 8000; height = 2000;
total_spots = 81;
total_anom = zeros(1, 81);
total_xvec = zeros(1, 81);

spots = 20; from = 50000; to = 2500; interval = -2500;
xvec = from:interval:to;
anom = zeros(1, spots);
for no = 1:spots
    rad = from + (no - 1)*interval;
    gap = 800;
    N = subheight/gap;
    anomaly = 0;
    for n = 0:N-1
       Radius = 2000 - (gap/4)*n;
       fc = subs(f, [z, x, h], [depth + gap*n, rad, 0]);
       func = matlabFunction(fc, 'Vars', {r, a});
       layer = integral2(func, 0, Radius, 0, 2*pi);
       anomaly = anomaly + dense * G * (depth + n*gap) * gap * layer;
    end
    anom(no) = anomaly;
end
anom = 10^5 * anom;
total_anom(1, 1:spots) = anom;
total_anom(1, (total_spots - spots + 1):total_spots) = fliplr(anom);
current = spots + 1;
total_xvec(1, 1:spots) = -xvec;
total_xvec(1, (total_spots - spots + 1):total_spots) = fliplr(xvec);

spots = 21; from = 2000; to = 0; interval = -100;
xvec = from:interval:to;
anom = zeros(1, spots);
for no = 1:spots
    rad = from + (no - 1)*interval;
    gap = 800;
    N = subheight/gap;
    anomaly = 0;
    elevation = (no - 1) * 2000 / (spots - 1);
    for n = 0:N-1
       Radius = 2000 - (gap/4)*n;
       fc = subs(f, [z, x, h], [depth + gap*n, rad, elevation]);
       func = matlabFunction(fc, 'Vars', {r, a});
       layer = integral2(func, 0, Radius, 0, 2*pi);
       anomaly = anomaly + dense * G * (depth + n*gap + elevation) * gap * layer;
    end
    anom(no) = anomaly;
end
anom = 10^5 * anom;
total_anom(1, current:(current + spots - 1)) = anom;
total_anom(1, (total_spots - current - spots + 2):(total_spots - current + 1)) = fliplr(anom);
total_xvec(1, current:(current + spots - 1)) = -xvec;
total_xvec(1, (total_spots - current - spots + 2):(total_spots - current + 1)) = fliplr(xvec);

subplot(2, 2, 1);
plot(total_xvec, total_anom);

  

% 生成地表物体引起的重力异常
% 为生成自由空气异常,需先执行计算布格重力异常的脚本(前)加以叠加
syms r a z x h;
f = r / sqrt((z - h)^2 + r^2 + x^2 - 2*r*x*cos(a))^3;

dense = 2800; G = 6.67E-11;
height = 2000;
total_spots = 81;
total_anom1 = zeros(1, 81);
total_xvec1 = zeros(1, 81);

spots = 20; from = 50000; to = 2500; interval = -2500;
xvec = from:interval:to;
anom = zeros(1, spots);
for no = 1:spots
    rad = from + (no - 1)*interval;
    gap = 200;
    N = height/gap;
    anomaly = 0;
    for n = 0:N-1
       Radius = 2000 - gap * (n + 0.5);
       fc = subs(f, [z, x, h], [gap*n + 100, rad, 0]);
       func = matlabFunction(fc, 'Vars', {r, a});
       layer = integral2(func, 0, Radius, 0, 2*pi);
       anomaly = anomaly - dense * G * n*gap * gap * layer;
    end
    anom(no) = anomaly;
end
anom = 10^5 * anom;
total_anom1(1, 1:spots) = anom;
total_anom1(1, (total_spots - spots + 1):total_spots) = fliplr(anom);
current = spots + 1;
total_xvec1(1, 1:spots) = -xvec;
total_xvec1(1, (total_spots - spots + 1):total_spots) = fliplr(xvec);

spots = 21; from = 2000; to = 0; interval = -100;
xvec = from:interval:to;
anom = zeros(1, spots);
for no = 1:spots
    rad = from + (no - 1)*interval;
    gap = 50;
    N = height/gap;
    anomaly = 0;
    elevation = (no - 1) * 2000 / (spots - 1);
    for n = 0:N-1
       Radius = 2000 - gap*(n + 0.5);
       fc = subs(f, [z, x, h], [gap*n + 25, rad, elevation + 2]);
       func = matlabFunction(fc, 'Vars', {r, a});
       layer = integral2(func, 0, Radius, 0, 2*pi);           
       anomaly = anomaly + dense * G * (elevation - n*gap - 25) * gap * layer;
    end
    anom(no) = anomaly;
end
anom = 10^5 * anom;
total_anom1(1, current:(current + spots - 1)) = anom;
total_anom1(1, (total_spots - current - spots + 2):(total_spots - current + 1)) = fliplr(anom);
total_xvec1(1, current:(current + spots - 1)) = -xvec;
total_xvec1(1, (total_spots - current - spots + 2):(total_spots - current + 1)) = fliplr(xvec);

subplot(2, 2, 2);
plot(total_xvec1, total_anom1);

freeair_xvec = total_xvec;
freeair = total_anom + total_anom1;
subplot(2, 2, 3); plot(freeair_xvec, freeair);

  

% 生成总重力异常
% 需要先执行布格重力异常脚本和自由空气异常脚本
total_spots = 81;
total_anom2 = zeros(1, 81);
total_xvec2 = zeros(1, 81);

spots = 20; from = 50000; to = 2500; interval = -2500;
xvec = from:interval:to;
total_xvec2(1, 1:spots) = -xvec;
total_xvec2(1, (total_spots - spots + 1):total_spots) = fliplr(xvec);

current = 21;
spots = 21; from = 2000; to = 0; interval = -100;
xvec = from:interval:to;
anom = zeros(1, spots);
for no = 1:spots
    elevation = (no - 1) * 2000 / (spots - 1);
    anom(no) = - elevation * 0.308;
end
total_anom2(1, current:(current + spots - 1)) = anom;
total_anom2(1, (total_spots - current - spots + 2):(total_spots - current + 1)) = fliplr(anom);
total_xvec2(1, current:(current + spots - 1)) = -xvec;
total_xvec2(1, (total_spots - current - spots + 2):(total_spots - current + 1)) = fliplr(xvec);

gravity_anomaly = freeair + total_anom2;

subplot(2, 2, 4); plot(freeair_xvec, gravity_anomaly);

  

转载于:https://www.cnblogs.com/gentle-min-601/p/10036685.html

  • 0
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值