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


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

    地壳密度 $kgcdot m^{-3}$ 上地幔密度 $gcdot cm^{-3}$ 地表地形圆锥体半径 (km) 地表地形圆锥体高度 (km) 计算莫霍面形变圆锥半径 (km) 计算莫霍面形变圆锥高度 (km) 地壳厚度 (km)
    $2.8 imes 10^{3}$ $3.5 imes 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);
    

      

  • 相关阅读:
    替换TStrings
    WordPress数据备份方案
    图像反色
    通过网络复制文件
    SQL Server的patindex和charindex的用法
    C冒泡排序 @100到200素数
    正则。IP验证
    C以二进制读、写、文本
    HTML下拉框、文本框、复选框!
    HTM页面获得本机时间
  • 原文地址:https://www.cnblogs.com/gentle-min-601/p/10036685.html
Copyright © 2020-2023  润新知