如何绘制以X、Y为坐标、z为时间的孔隙度3D变异函数？

下面针对你使用过的GSTools、GStat工具，给出具体的3D变异函数（X/Y为空间坐标、Z为时间维度）实现步骤和代码：

一、GSTools实现（推荐，对多维变异函数支持更灵活） #

GSTools原生支持3维变异函数计算，可直接将时间维度作为第三维处理。

1. 环境与数据准备 #

先导入依赖库，替换模拟数据为你的真实孔隙度数据集：

import numpy as np
import gstools as gs
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 替换为你的真实数据：X/Y坐标、时间Z、孔隙度porosity
np.random.seed(42)
x = np.random.rand(1000) * 100  # X空间坐标（示例范围0-100）
y = np.random.rand(1000) * 100  # Y空间坐标（示例范围0-100）
z = np.random.rand(1000) * 50   # 时间维度（示例范围0-50）
porosity = np.random.randn(1000) * 0.05 + 0.3  # 孔隙度数据

2. 计算实验变异函数 #

初始化3维变异函数对象，设置滞后距区间（bins）：

# 定义滞后距区间，可根据数据尺度调整
bins = np.linspace(0, 100, 30)
# 计算3D实验变异函数
vario = gs.Variogram(
    pos=[x, y, z],  # 三维位置：X/Y/时间Z
    field=porosity,
    bins=bins,
    dim=3,
    normalizer=None  # 若数据分布偏态，可设置归一化器如gs.normalizer.BoxCox
)

3. 拟合理论变异函数模型 #

支持球状、指数、高斯等多种模型，这里以球状模型为例：

# 拟合3维球状模型
vario.fit_model(model=gs.Spherical(dim=3))
# 打印拟合得到的参数（变程、基台值、块金值）
print("拟合模型参数：", vario.model)

4. 可视化 #

3D变异函数云图（展示所有数据对的滞后距与半方差）

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 提取滞后距分量与对应半方差
lag_x, lag_y, lag_z, gamma = vario.cloud
scatter = ax.scatter(lag_x, lag_y, lag_z, c=gamma, cmap='viridis', alpha=0.6)
ax.set_xlabel('X方向滞后距')
ax.set_ylabel('Y方向滞后距')
ax.set_zlabel('时间滞后距')
ax.set_title('孔隙度3D变异函数云图')
plt.colorbar(scatter, label='半方差')
plt.show()

径向变异函数（各向同性假设下，滞后距长度与半方差的关系）

fig, ax = plt.subplots(figsize=(8, 6))
# 实验变异函数点
ax.plot(vario.bins, vario.gamma, 'o', label='实验变异函数')
# 拟合的理论模型曲线
ax.plot(vario.bins, vario.model(vario.bins), '-', label='拟合球状模型')
ax.set_xlabel('滞后距长度')
ax.set_ylabel('半方差')
ax.set_title('孔隙度径向变异函数（空间+时间联合）')
ax.legend()
plt.grid(True)
plt.show()

5. 各向异性设置（可选） #

若空间与时间维度的变异特性不同（如时间变程远小于空间变程），可设置各向异性参数：

# 定义各向异性模型：X/Y方向变程50，时间方向变程25
model = gs.Spherical(
    dim=3, 
    len_scale=[50, 50, 25],  # 三个维度的变程
    anisotropy=[1, 0.5],     # 各向异性比值
    angles=[0, 0, 0]         # 旋转角度（无旋转则设为0）
)
vario.fit_model(model=model)

二、GStat实现 #

GStat同样支持3维变异函数计算，步骤如下：

import gstat as gs_g
import matplotlib.pyplot as plt

# 计算3D实验变异函数
vario_g = gs_g.Variogram(
    coordinates=np.vstack((x, y, z)).T,
    values=porosity,
    n_lags=30,
    maxlag=100,
    dim=3
)

# 拟合球状模型
vario_g.fit(model='spherical')

# 可视化实验变异函数与拟合模型
vario_g.plot(show=False)
plt.title('GStat 3D孔隙度变异函数')
plt.show()

关键注意事项 #

• 尺度归一化：若时间与空间的单位/尺度差异较大，建议先对时间维度做归一化（如缩放至与空间尺度相当），避免变异函数受尺度偏差影响。
• 数据量要求：3D变异函数需要足够的数据对支撑，否则实验变异函数会存在严重噪声，建议数据集规模不小于1000条。
• 模型选择：可通过拟合优度（如R²）选择最优理论模型，GSTools与GStat均支持模型评估。

内容的提问来源于stack exchange，提问作者Noor Ul huda Choudhary