4. 实验变异函数计算实战

实验变异函数是连接“样品数据”与“理论模型”的桥梁。它的计算并不复杂，但真正决定结果可用性的，是滞后距分组方式、距离容差与角度容差的设置，以及如何从云图和箱线图中读出趋势与漂移。本节以一个假想的斑岩铜矿钻孔品位数据集为例，完整走一遍从数据准备到趋势识别的流程。示例中的参数与数值均为教学假设，用于说明操作逻辑与判读依据，并非实测结果。

一、输入数据与坐标约定

假设我们有一个包含 3 个钻孔、共 480 个样的品位数据集，字段为：孔号 BH_ID、样号 SAMP_ID、三维坐标 X/Y/Z（X 向东、Y 向北、Z 向上，单位 m）、样长 LENGTH（多为 1.0 m，少量 0.5 m）、品位 CU（%）。计算前必须明确三点：坐标参照系与单位统一；样长不等时需按长度加权；对高度偏倚的品位（铜品位常近似对数正态）通常先做对数变换再计算变异函数，以避免少数高品位样主导结果。

数据清洗的常见动作包括：剔除 CU 为负或缺失的记录；检查坐标重复点；对超出地质合理上限的极值（如 >5%）单独标记，决定是否截尾（如截尾至 99% 分位）。这些处理会直接改变实验变异函数的形态，必须在报告中记录。

二、滞后距分组：把“点对”变成“距离类”

实验变异函数的本质是对所有样品点对按距离分箱，计算每箱内半方差均值。设滞后距步长 lag = 25 m，滞后距容差取步长的一半（lag_tol = 12.5 m），最大距离取勘探网度的 1/2～2/3，示例取 max_dist = 400 m。于是第 k 个滞后距中心为 25k，凡点对距离落在 [25k−12.5, 25k+12.5) 内即归入该箱。全向计算时角度容差设为 90°，即不区分方向，所有点对参与统计。

每个滞后距箱需同时记录：点对数 N(k)、平均距离 h(k)、半方差均值 γ(h)。点对数过少（经验上少于 30 对）的箱应合并或舍弃，否则该点会剧烈跳动，误导后续拟合。半方差按经典公式计算：

γ(h) = 1 / (2·N(h)) · Σ [ Z(xi) − Z(xi+h) ]²
其中 N(h) 为距离落入 h 箱的点对数，
Z 为品位（或其对数值），xi 与 xi+h 为点对两端坐标。

若样长不等，则用长度加权形式，以 wi·wj 作为点对权重，分母相应改为 2·Σwi·wj。这一步在软件中通常以“长度加权”选项体现，务必与数据实际情况匹配。

三、变异函数云图：先看“散”，再看“势”

变异函数云图（variogram cloud）把每一个点对的半方差 0.5·(Zi−Zj)² 对距离 h 作散点图，不做分箱平均。它的价值在于暴露被平均值掩盖的结构：如果云图在某一距离后整体抬升，说明存在明显的空间变异性增大；如果云图呈现上下两条“分支”，往往意味着数据中存在趋势或不同岩性域混合；若近距离出现大量高值散点，则提示可能存在短尺度高品位簇或数据错误。

云图点对数量随样品数呈平方增长（480 个样约 11.5 万对），直接绘制会非常密集。实操中常做两件事：一是对半方差取对数纵轴，压缩高值区；二是按距离抽样或按分位数着色，突出主体分布。云图只用于定性诊断，不直接用于拟合。

四、滞后距分组箱线图：定量看每个距离的分布

箱线图是对云图的分箱汇总：每个滞后距箱画出一个箱体，箱内横线为中位数，箱体上下边为四分位数，须线延伸至 1.5 倍四分位距内的极值，箱外点为离群点对。相比单一均值曲线，箱线图同时给出该距离上半方差的集中程度与离散程度。判读要点如下：

• 中位数随距离上升并趋于平稳：典型的变异性随距离增大后达到基台，说明存在可拟合的变程。
• 箱体随距离明显变宽：远距离点对少且混入不同区域，统计不稳，拟合时应降低这些箱的权重或截断最大距离。
• 近距离箱体中位数异常偏高：可能存在块金效应显著，或样品间距小于采样尺度导致重复样、混样问题。
• 各箱中位数整体呈单调倾斜而非水平：提示存在大尺度趋势（漂移），需先去除趋势再计算残差变异函数。

五、从图中识别趋势与漂移

趋势（trend）指品位随坐标呈系统性变化，漂移（drift）是趋势在变异函数上的表现：当存在漂移时，实验变异函数不会在有限距离内趋于平稳基台，而是随距离持续上升，甚至呈抛物线状。判读时把云图与箱线图对照看：若箱线图中位数曲线在最大滞后距处仍明显上升、未出现平台，且云图远距离散点整体抬升并出现分叉，就应怀疑存在大尺度趋势。

确认趋势的常用做法是：分别沿 X、Y、Z 方向计算方向性实验变异函数，若某方向在远距离显著高于其他方向且不收敛，说明该方向存在漂移；也可直接对品位与坐标做线性或二次回归，若回归系数显著且残差变异函数收敛，则用残差建模、再把趋势加回。示例中若发现 Z 方向（垂向）品位随深度系统升高，就应先拟合垂向趋势面，再对残差计算全向变异函数。

六、一份可直接套用的计算记录表

为了让计算过程可复现、可审查，建议每个滞后距箱都留一行记录。下表为示例结构（数值为教学假设）：

记录表的作用不只是存档：当后续拟合球状或指数模型时，可直接用 N(h) 作为加权依据，让点对多的近距离箱主导拟合，避免远距离稀疏箱把变程拉偏。若某箱 N(h) 低于 30，应在备注中标注“合并或剔除”，并在最终拟合时排除。

七、常见问题与判断依据

第一，块金值偏大：若 h=25 m 箱的中位数已接近基台的一半，说明短尺度变异性强，可能是采样间距大于矿化尺度，或存在分析误差，此时不应强行把块金压小。第二，曲线抖动：多因分组过细或点对不足，可增大 lag 或合并相邻箱。第三，远距离翘尾：优先排查趋势与不同岩性域混合，而不是直接加大基台。第四，对数变换后结果需还原：在对数空间拟合的变程可直接使用，但基台与块金在还原到原始品位尺度时需谨慎解释，报告中应注明计算所用变换。

完成上述步骤后，你得到的不是一条“好看的曲线”，而是一组带点对数、带分组参数、带趋势判断的实验变异函数结果。它将成为后续各向异性分析与理论模型拟合的输入，也是判断“数据是否支持地质统计建模”的第一道证据。