告别“彩虹噪音”：地球物理勘探中的三维密度散点图绘制之道
三维密度散点图？别闹了！
各位，别一上来就想着 scatter3 怎么用，先扪心自问一句：你真的理解“密度”的意义吗？现在网上那些教程，张口闭口就是“三维密度散点图”，代码复制粘贴一套一套的，画出来的图呢？一堆彩色的点，美其名曰“密度分布”，实际上就是毫无物理意义的噪音！
地球物理数据，那可不是随便生成的随机数。它背后代表的是地下的秘密：矿物富集程度、岩石破碎带的分布、污染物扩散浓度……没有这些概念，你画“密度”给谁看？给外星人吗？
数据准备：巧妇难为无米之炊
地球物理数据的特殊性，决定了数据准备的重要性。想象一下，你拿到了一份矿产勘探数据，里面包含经纬度、深度，以及对应的品位值（比如金的含量，单位：g/t）。要画三维密度散点图，第一步是什么？当然不是直接丢进 scatter3 里！
- 数据格式转换： 将数据整理成MATLAB可以处理的矩阵形式。例如，创建一个N×4的矩阵，其中N是数据点的数量，四列分别是经度、纬度、深度和品位值。
- 坐标转换： 如果你的数据是经纬度坐标，可能需要将其转换为投影坐标系（例如UTM），以避免在计算距离和面积时出现偏差。MATLAB有相应的工具箱可以进行坐标转换。
- 数据插值： 原始数据点通常是离散的，为了获得更平滑的密度分布，可能需要进行数据插值。常用的插值方法包括线性插值、立方插值等。但要注意，插值过度可能会掩盖真实的数据特征。
案例： 假设我们在某金矿区采集了钻孔岩心数据，获得了每个岩心样品的金品位和空间位置。为了分析金的富集规律，我们需要将这些数据整理成MATLAB可以处理的格式，并进行适当的插值。
密度计算：别再用那些“想当然”的函数了！
kde2d？histcounts？这些函数用起来是方便，但你考虑过它们的局限性吗？地球是球体，不是平的！这些函数考虑了地球曲率吗？地球物理数据通常是非均匀分布的，它们能处理吗？
是时候祭出一些更高级的武器了！
反距离权重法 (IDW)：简单实用，居家旅行必备
反距离权重法 (Inverse Distance Weighting, IDW) 的原理很简单：距离越近，权重越大。对于某个数据点，它的密度可以定义为周围数据点品位的加权平均值，权重与距离的负幂次成正比。
公式：
$D_i = \frac{\sum_{j=1}^{N} \frac{w_j}{d_{ij}^p}}{\sum_{j=1}^{N} \frac{1}{d_{ij}^p}}$
其中，$D_i$ 是点 i 的密度估计值，$w_j$ 是点 j 的品位值，$d_{ij}$ 是点 i 和点 j 之间的距离，$p$ 是距离衰减指数 (通常取1或2)。
MATLAB代码示例：
function density = idw_density(x, y, z, values, xi, yi, zi, p)
% IDW_DENSITY 使用反距离权重法估计三维散点的密度
% x, y, z: 原始数据点的坐标
% values: 原始数据点的品位值
% xi, yi, zi: 需要估计密度的点的坐标
% p: 距离衰减指数
n = length(x);
sum_weights = 0;
weighted_sum = 0;
for j = 1:n
distance = sqrt((xi - x(j))^2 + (yi - y(j))^2 + (zi - z(j))^2);
if distance == 0
density = values(j); % 如果距离为0，直接取该点的值
return;
end
weight = 1 / (distance^p);
sum_weights = sum_weights + weight;
weighted_sum = weighted_sum + weight * values(j);
end
density = weighted_sum / sum_weights;
end
克里金法 (Kriging)：效果杠杠的，就是有点复杂
克里金法是一种更高级的地统计学方法，它不仅考虑了距离，还考虑了数据的空间自相关性。克里金法通过构建变异函数来描述数据的空间相关性，并利用变异函数来估计密度。虽然克里金法计算量更大，但通常能提供更准确的密度估计。
简单来说，克里金法就是： 先研究数据之间的“亲疏关系”（空间自相关性），再根据这种关系来进行插值。就像相亲一样，先了解对方的家庭背景、性格爱好，再决定是否发展关系。
MATLAB代码示例 (简化版)：
% 注意：这只是一个简化的示例，实际应用中需要使用更完善的克里金工具箱
function density = simple_kriging(x, y, z, values, xi, yi, zi, variogram_model, range, sill)
% SIMPLE_KRIGING 使用简单克里金法估计三维散点的密度
% x, y, z: 原始数据点的坐标
% values: 原始数据点的品位值
% xi, yi, zi: 需要估计密度的点的坐标
% variogram_model: 变异函数模型 (例如 'spherical', 'exponential')
% range: 变异函数的作用范围
% sill: 变异函数的基台值
n = length(x);
A = zeros(n + 1, n + 1);
b = zeros(n + 1, 1);
% 构建克里金矩阵
for i = 1:n
for j = 1:n
distance = sqrt((x(i) - x(j))^2 + (y(i) - y(j))^2 + (z(i) - z(j))^2);
A(i, j) = variogram(distance, variogram_model, range, sill); % 使用变异函数计算协方差
end
A(i, n + 1) = 1;
A(n + 1, i) = 1;
end
A(n + 1, n + 1) = 0;
% 构建右侧向量
for i = 1:n
distance = sqrt((xi - x(i))^2 + (yi - y(i))^2 + (zi - z(i))^2);
b(i) = variogram(distance, variogram_model, range, sill);
end
b(n + 1) = 1;
% 求解克里金方程
weights = A \ b;
% 计算密度估计值
density = sum(weights(1:n) .* values);
end
function gamma = variogram(h, model, range, sill)
% VARIOGRAM 计算变异函数值
switch lower(model)
case 'spherical'
if h <= range
gamma = sill * (1.5 * (h / range) - 0.5 * (h / range)^3);
else
gamma = sill;
end
case 'exponential'
gamma = sill * (1 - exp(-3 * h / range));
otherwise
error('Unsupported variogram model.');
end
end
自适应带宽核密度估计：灵活应变，效果更佳
固定带宽的核密度估计 (Kernel Density Estimation, KDE) 就像是用一把固定大小的刷子涂颜色，对于密度变化大的区域，可能会出现过度平滑或欠平滑的问题。自适应带宽KDE则可以根据数据密度自动调整带宽，在密度高的区域使用较小的带宽，在密度低的区域使用较大的带宽，从而获得更精细的密度估计结果。
三维散点图绘制：告别“彩虹色”，拥抱科学之美
颜色映射：别再用“彩虹色”了，审美呢？
默认的“彩虹色”颜色映射（jet）对人眼极不友好，容易产生视觉偏差。相邻颜色之间的差异过大，会让人误以为数据之间存在明显的跳跃。推荐使用更科学的颜色映射方案，例如 parula、viridis 或者自定义的颜色映射。
colormap(parula); % 使用parula颜色映射
% 或者
colormap(viridis); % 使用viridis颜色映射
颜色标尺 (Colorbar)：让数据说话
颜色标尺是三维散点图的灵魂。没有颜色标尺，谁知道颜色代表什么？要设置颜色标尺的范围和标签，使其能够清晰地反映数据的物理意义。例如，如果颜色代表金品位，可以将颜色标尺的范围设置为0到最大品位值，并在颜色标尺上标注品位值的单位（g/t）。
c = colorbar;
c.Label.String = 'Gold Grade (g/t)';
caxis([0, max(values)]); % 设置颜色范围
透明度 (Alpha)：让密度“浮出水面”
使用透明度可以增强三维散点图的可视化效果。可以根据密度值调整点的透明度，使高密度区域更加突出。例如，可以将密度值较高的点的透明度设置为较低，使其更加显眼。
alpha = density / max(density); % 将密度值归一化到0-1之间
scatter3(x, y, z, 50, values, 'filled', 'MarkerFaceAlpha', 0.5); % 设置透明度
视角调整：找到最佳观察角度
使用 view 函数调整视角，找到最佳的观察角度。不同的视角可以揭示数据的不同特征。例如，可以从正上方观察数据的平面分布，也可以从侧面观察数据的垂直分布。
view([45 30]); % 设置视角
完整代码示例
% 模拟地球物理勘探数据
num_points = 1000;
x = rand(num_points, 1) * 1000; % x坐标
y = rand(num_points, 1) * 1000; % y坐标
z = rand(num_points, 1) * 100; % 深度坐标
values = 5 * exp(-((x - 500).^2 + (y - 500).^2 + (z - 50).^2) / (2 * 200^2)) + randn(num_points, 1); % 模拟品位值
% 密度估计 (使用IDW)
xi = linspace(min(x), max(x), 50);
yi = linspace(min(y), max(y), 50);
zi = linspace(min(z), max(z), 20);
[X, Y, Z] = meshgrid(xi, yi, zi);
density = zeros(size(X));
for i = 1:numel(X)
density(i) = idw_density(x, y, z, values, X(i), Y(i), Z(i), 2); % 使用IDW计算密度
end
% 三维散点图绘制
figure;
scatter3(X(:), Y(:), Z(:), 50, density(:), 'filled');
% 设置颜色映射
colormap(parula);
% 设置颜色标尺
c = colorbar;
c.Label.String = 'Simulated Value';
caxis([min(density(:)), max(density(:))]);
% 设置坐标轴标签和标题
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Depth (m)');
title('3D Density Scatter Plot of Simulated Geophysical Data');
% 调整视角
view([45 30]);
% 添加透明度
alpha = density(:) / max(density(:));
scatter3(X(:), Y(:), Z(:), 50, density(:), 'filled', 'MarkerFaceAlpha', 0.5);
%IDW函数代码
function density = idw_density(x, y, z, values, xi, yi, zi, p)
% IDW_DENSITY 使用反距离权重法估计三维散点的密度
% x, y, z: 原始数据点的坐标
% values: 原始数据点的品位值
% xi, yi, zi: 需要估计密度的点的坐标
% p: 距离衰减指数
n = length(x);
sum_weights = 0;
weighted_sum = 0;
for j = 1:n
distance = sqrt((xi - x(j))^2 + (yi - y(j))^2 + (zi - z(j))^2);
if distance == 0
density = values(j); % 如果距离为0，直接取该点的值
return;
end
weight = 1 / (distance^p);
sum_weights = sum_weights + weight;
weighted_sum = weighted_sum + weight * values(j);
end
density = weighted_sum / sum_weights;
end
案例分析：矿产勘探中的应用
假设我们使用上述代码处理了某铜矿区的钻孔数据，得到了一个三维密度散点图。通过观察该图，我们可以发现高密度区域主要集中在地下100米至200米之间，且呈现出明显的条带状分布。这可能表明该矿区存在一个埋藏较深的层状矿体，且矿体的走向与条带状分布的方向一致。结合地质资料，我们可以进一步推断矿体的成因和规模，为后续的勘探工作提供指导。
总结与展望
三维密度散点图是地球物理数据可视化的有力工具。但要绘制出真正有价值的密度散点图，需要深入理解数据的物理意义，选择合适的密度估计方法，并掌握科学的颜色映射和视角调整技巧。未来，随着虚拟现实 (VR) 和增强现实 (AR) 等技术的不断发展，三维密度散点图的可视化效果将得到进一步提升，为地球物理勘探带来更广阔的应用前景。
别再满足于那些“千篇一律”的教程了！用你的智慧和创造力，去探索地球深处的奥秘吧！2026年，地球物理的未来，掌握在你们手中！