用MATLAB解决文件中的问题

要使用MATLAB完成《勘探地球物理学实验指导书》中的实验任务，可以按照以下步骤编写代码：

### 实验一：反射地震波合成与动校正处理

#### 1. 计算反射系数和传播时间

% 参数设置
rho1 = 2500; % 上层介质密度 (kg/m^3)
rho2 = 3000; % 下层介质密度 (kg/m^3)
v1 = 3000;   % 上层介质速度 (m/s)
v2 = 4000;   % 下层介质速度 (m/s)
d = 1000;    % 层厚度 (m)
xi = 0:25:1000; % 炮检距 (m)

% 计算反射系数
R = (rho2 * v2 - rho1 * v1) / (rho2 * v2 + rho1 * v1);

% 计算反射时间
t0 = 2 * d / v2;
ti = sqrt(t0^2 + (xi.^2) / (v2^2));

#### 2. 合成地震道

% 子波参数
fm = 30; % 主频 (Hz)
dt = 0.002; % 时间采样间隔 (s)
t = -0.5:dt:0.5; % 时间轴 (s)
w = (1 - 2 * pi^2 * fm^2 * t.^2) .* exp(-pi^2 * fm^2 * t.^2); % Ricker子波

% 合成地震道
n = length(ti);
seis = zeros(n, length(t));
for i = 1:n
    seis(i, :) = conv(R * dirac(t - ti(i)), w, 'same');
end

#### 3. 动校正处理

% 正常时差
NMO = @(xi, v, t0) sqrt(t0^2 + (xi.^2) / (v^2)) - t0;

% 动校正
seis_nmo = zeros(size(seis));
for i = 1:n
    nmo_i = NMO(xi(i), v2, t0);
    idx = round(nmo_i / dt) + 1;
    if idx <= size(seis, 2)
        seis_nmo(i, 1:end-idx+1) = seis(i, idx:end);
    end
end

% 显示结果
figure;
imagesc(t, xi, seis);
colorbar;
title('原始地震道');
xlabel('时间 (s)');
ylabel('炮检距 (m)');

figure;
imagesc(t, xi, seis_nmo);
colorbar;
title('动校正后的地震道');
xlabel('时间 (s)');
ylabel('炮检距 (m)');

### 实验二：重磁异常正演与异常分离

#### 1. 计算球体重力异常

% 参数设置
G = 6.67e-11; % 万有引力常数 (m^3/kg/s^2)
R1 = 100;     % 球体1半径 (m)
R2 = 150;     % 球体2半径 (m)
rho1 = 500;   % 球体1剩余密度 (kg/m^3)
rho2 = 700;   % 球体2剩余密度 (kg/m^3)
h1 = 500;     % 球体1埋深 (m)
h2 = 1000;    % 球体2埋深 (m)
x1 = 500;     % 球体1中心 x 坐标 (m)
x2 = 900;     % 球体2中心 x 坐标 (m)
x = -1000:10:1000; % 测点 x 坐标 (m)

% 重力异常
g1 = G * rho1 * R1^3 * (x - x1) ./ ((x - x1).^2 + h1^2).^(3/2);
g2 = G * rho2 * R2^3 * (x - x2) ./ ((x - x2).^2 + h2^2).^(3/2);
g_total = g1 + g2;

% 显示结果
figure;
plot(x, g_total);
hold on;
plot(x, g1, '--');
plot(x, g2, ':');
legend('总重力异常', '球体1重力异常', '球体2重力异常');
xlabel('测点 x 坐标 (m)');
ylabel('重力异常 (mGal)');
title('重力异常曲线');

#### 2. 计算球体磁异常

% 参数设置
mu0 = 4 * pi * 1e-7; % 真空磁导率 (H/m)
M1 = 1000;           % 球体1磁化强度 (Am^2)
M2 = 1500;           % 球体2磁化强度 (Am^2)
V1 = 4/3 * pi * R1^3; % 球体1体积 (m^3)
V2 = 4/3 * pi * R2^3; % 球体2体积 (m^3)
m1 = M1 * V1;         % 球体1磁矩 (Am^2)
m2 = M2 * V2;         % 球体2磁矩 (Am^2)
I1 = 60;              % 球体1磁倾角 (度)
I2 = 70;              % 球体2磁倾角 (度)
A1 = 0;               % 球体1磁偏角 (度)
A2 = 0;               % 球体2磁偏角 (度)

% 转换角度为弧度
I1_rad = deg2rad(I1);
I2_rad = deg2rad(I2);
A1_rad = deg2rad(A1);
A2_rad = deg2rad(A2);

% 磁异常
B1 = mu0 / (4 * pi) * m1 * (cos(I1_rad) * (x - x1) + sin(I1_rad) * cos(A1_rad) * h1) ./ ...
      ((x - x1).^2 + h1^2).^(3/2);
B2 = mu0 / (4 * pi) * m2 * (cos(I2_rad) * (x - x2) + sin(I2_rad) * cos(A2_rad) * h2) ./ ...
      ((x - x2).^2 + h2^2).^(3/2);
B_total = B1 + B2;

% 显示结果
figure;
plot(x, B_total);
hold on;
plot(x, B1, '--');
plot(x, B2, ':');
legend('总磁异常', '球体1磁异常', '球体2磁异常');
xlabel('测点 x 坐标 (m)');
ylabel('磁异常 (nT)');
title('磁异常曲线');

#### 3. 异常分离（向上延拓）

% 延拓高度
dh = 100; % 延拓高度 (m)

% 延拓系数
C0 = 0.2952;
C1 = 0.1652;
C2 = 0.0660;
C3 = 0.0325;
C4 = 0.0190;
C5 = 0.0124;
C6 = 0.0087;

% 上延拓
g_up = C0 * g_total + C1 * circshift(g_total, [0, 1]) + C2 * circshift(g_total, [0, 2]) + ...
       C3 * circshift(g_total, [0, 3]) + C4 * circshift(g_total, [0, 4]) + ...
       C5 * circshift(g_total, [0, 5]) + C6 * circshift(g_total, [0, 6]);

% 显示结果
figure;
plot(x, g_total, '-');
hold on;
plot(x, g_up, '--');
legend('原重力异常', '上延拓后重力异常');
xlabel('测点 x 坐标 (m)');
ylabel('重力异常 (mGal)');
title('重力异常上延拓');

以上代码涵盖了《勘探地球物理学实验指导书》中的主要实验内容，并提供了详细的实现步骤。你可以根据具体需求调整参数和代码。希望这些代码对你有所帮助！