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

实验变异函数是普通克里金插值中连接“样品数据”与“理论模型”的桥梁。它不依赖任何先验假设，而是直接从组合品位样品的空间配对中统计出品位随距离变化的实际规律。本节以某铁矿床的 186 个组合品位样品为例（示例数据，用于教学演示），完整走一遍从样品对生成、滞后距分组、角度容差筛选到多方向实验变异函数曲线输出的过程，并说明每一步参数选择对最终曲线形态的影响。

一、计算对象与基本公式

组合品位样品是把相邻工程或相邻品位的原始样品按长度加权合并后的结果，其空间支撑更大、品位更平滑，能有效抑制特高品位对变异函数的干扰。设样品位置为 u_i，组合品位为 Z(u_i)，则实验变异函数的基本估计量为：

γ*(h) = 1 / (2N(h)) · Σ [ Z(u_i) − Z(u_i + h) ]²

其中 h 为滞后距向量，N(h) 为落在该滞后距容差范围内的样品对数量。注意公式中除以 2N(h)，是因为每一对样品 (i, j) 与 (j, i) 被重复计入，取半即为变异函数的无偏估计。实际编程时通常只遍历 i < j 的样品对，此时分母直接用 N(h) 即可，两种写法等价，但必须前后统一，否则曲线整体会被放大或缩小一倍。

二、滞后距与容差的确定

滞后距 h 的步长和最大距离决定了曲线的分辨率与可信区间。经验做法是：最大滞后距取研究区最大样品间距的一半左右，本例样品分布范围约 1200 m，故取 h_max = 600 m；步长取最小平均样品间距的量级，本例约 50 m，于是得到 12 个滞后距级别：50、100、150……600 m。步长过小会导致每级样品对过少、曲线剧烈跳动；步长过大则会把短距离的结构性变化抹平。

容差 Δh 用于把“恰好等于 h”的严格条件放宽为一个区间 [h − Δh, h + Δh]。常用取法是 Δh = 步长 / 2，本例即 ±25 m，这样相邻滞后距区间首尾相接、不重叠也不留空。若样品稀少，可放宽到步长的 0.6～0.75 倍以增加样品对数量，但需在报告中注明，因为它会轻微平滑曲线。

三、角度容差与方向划分

各向异性矿床的变异函数在不同方向上形态差异明显，因此需要按方向分别统计。方向由方位角 θ 定义，以正北为 0°、顺时针为正（本课程全文统一采用此约定）。角度容差 Δθ 表示样品对连线方向与目标方向夹角在 ±Δθ 内即被采纳。本例设置四个方向：0°（南北）、45°（北东）、90°（东西）、135°（北西），角度容差取 22.5°，恰好把 180° 平面均分为四个互不重叠的扇区。

需要特别提醒：变异函数具有对称性，γ(h) = γ(−h)，因此方向只需在 0°～180° 内划分，0° 与 180° 是同一方向。若角度容差设置过大（如 45°），扇区之间会重叠，同一对样品被计入多个方向，导致曲线被人为拉平；设置过小（如 10°）则每个方向样品对不足，曲线噪声大。22.5° 是四方向划分下的常用折中值。

四、计算过程与结果判读

按上述参数执行后，每个方向、每个滞后距级别都会得到一个 N(h) 和一个 γ*(h)。以 0° 方向第一级 h = 50 m 为例，若落入 [25, 75] m 且方位角在 −22.5°～22.5° 内的样品对共 34 对，各对品位差的平方和为 41.2（品位单位为 %TFe，平方后单位为 %²），则 γ*(50) = 41.2 / 34 ≈ 1.21 %²。逐级计算即可得到该方向的离散点列。

把四个方向的点列画在同一坐标系中，就得到本节的核心输出——多方向实验变异函数曲线。判读时重点关注三点：其一，曲线是否随 h 增大而上升并趋于平稳，若持续上扬说明 h_max 取得过大或存在趋势，需先做趋势剔除；其二，各方向曲线达到平稳的滞后距（变程）是否一致，本例东西向约 450 m、南北向约 250 m，说明存在明显的几何各向异性，后续理论模型拟合应采用各向异性椭球；其三，短距离段曲线是否出现异常高值，这往往提示存在 nugget 效应或特高品位残留，应回查组合品位数据。

五、可直接套用的参数记录表

为保证计算可复现，建议在每次计算时同步填写下表（本例为示例填写值），它既是计算日志，也是后续理论模型拟合的输入依据。

完成本节计算后，得到的四组 (h, γ*(h), N(h)) 数据就是下一步理论变异函数模型拟合的直接输入。参数选择没有唯一正确答案，但每一次调整都应在记录表中留痕，这样当拟合残差偏大时，才能快速回溯是滞后距、容差还是方向划分造成的影响。