在完全均衡的模型下，若地表有一圆锥体（山峰等），计算跨越山顶的截面上所得到的各种重力异常。
| 地壳密度 $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;endanom = 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;endanom = 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;endanom = 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;endanom = 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;endtotal_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);