土壤属性插值复现

Author

土壤属性数字制图复现小组

Published

May 28, 2026

研究背景

本次作业复现普通克里格（Ordinary Kriging, OK） 空间插值方法，基于湖北省随州市曾都区654个表层土壤样点数据，实现土壤属性的连续空间分布制图，验证经典地统计方法在县域土壤属性数字制图中的应用效果。

数据预处理

数据源说明

土壤样点数据：曾都历史数据.xls，共654个样点，包含X/Y投影坐标、土壤属性（碱解氮、有效磷、速效钾、有机质、PH值）及地形、气候、土壤类型等环境变量
插值网格数据：曾都渔网点.shp，研究区规则网格点，用于生成连续空间分布结果

数据读取与预处理

# ====================== 创建processed文件夹 ======================
os.makedirs(PROCESSED_DIR, exist_ok=True)
print(f"✅ 已创建处理后数据文件夹：{PROCESSED_DIR}")

# ====================== 读取并预处理土壤样点数据 ======================
print("\n📥 正在读取原始土壤样点数据...")
try:
    # 读取Excel数据
    raw_df = pd.read_excel(RAW_SAMPLE_PATH, sheet_name="曾都历史数据")
    print(f"📊 原始样点数据：")
    print(f"   - 总记录数：{len(raw_df)}")
    print(f"   - 核心列名：{[X_COL, Y_COL, TARGET_ATTR, 'PH值', '碱解氮', '有效磷', '速效钾']}")
except Exception as e:
    print(f"❌ 读取Excel数据失败：{e}")
    exit(1)

# 3.1 数据清洗：删除缺失值
print("\n🧹 数据清洗：删除缺失值...")
processed_df = raw_df.dropna(subset=[X_COL, Y_COL, TARGET_ATTR])
print(f"   - 清洗后记录数：{len(processed_df)}")

# 3.2 去重
print("🗑️ 去除重复数据...")
initial_count = len(processed_df)
processed_df = processed_df.drop_duplicates(subset=[X_COL, Y_COL, TARGET_ATTR], keep="first")
print(f"   - 去除重复数：{initial_count - len(processed_df)}，剩余记录数：{len(processed_df)}")

# 3.3 异常值处理
print("🚫 处理异常值...")
values = processed_df[TARGET_ATTR]
mean_val = values.mean()
std_val = values.std()
lower_bound = mean_val - 3 * std_val
upper_bound = mean_val + 3 * std_val

initial_count = len(processed_df)
processed_df = processed_df[(processed_df[TARGET_ATTR] >= lower_bound) & 
                             (processed_df[TARGET_ATTR] <= upper_bound)]
print(f"   - 删除异常值：{initial_count - len(processed_df)} 条，剩余记录数：{len(processed_df)}")

# 3.4 转换为GeoDataFrame
print("\n🌐 转换为空间数据格式...")
gdf = gpd.GeoDataFrame(
    processed_df,
    geometry=gpd.points_from_xy(processed_df[X_COL], processed_df[Y_COL]),
    crs="EPSG:4547"  
)

gdf[f"pred_{TARGET_ATTR}"] = None
gdf.rename(columns={TARGET_ATTR: "organic"}, inplace=True) 

# ====================== 读取并预处理渔网点数据 ======================
print("\n📥 正在读取插值网格数据...")
try:
    grid_gdf = gpd.read_file(RAW_GRID_PATH)
    print(f"📊 网格数据信息：")
    print(f"   - 总网格点数：{len(grid_gdf)}")
    print(f"   - 坐标系：{grid_gdf.crs}")
except Exception as e:
    print(f"❌ 读取网格数据失败：{e}")
    exit(1)

# 4.1 确保网格和样点坐标系一致
if grid_gdf.crs != gdf.crs:
    print("🔄 统一网格坐标系与样点数据...")
    grid_gdf = grid_gdf.to_crs(gdf.crs)

普通克里格插值

普通克里格是地统计中最经典的无偏最优插值方法，假设区域化变量满足二阶平稳性，通过半方差函数拟合样点间的空间自相关性，实现未知点的属性预测。

# ===================== 主程序 =====================
if __name__ == '__main__':
    print("===== 开始克里格插值计算 =====")
    
    # 1. 读取数据
    sample_df = pd.read_csv(SAMPLE_PATH, encoding='utf-8')
    grid_gdf = gpd.read_file(GRID_PATH, encoding='gbk')

    # 2. 克里格插值
    x = sample_df["X"].values
    y = sample_df["Y"].values
    z = sample_df[TARGET_ATTR].values

    ok_model = OrdinaryKriging(
        x, y, z,
        variogram_model="gaussian",
        verbose=False,
        enable_plotting=False
    )

    grid_x = grid_gdf["X"].values
    grid_y = grid_gdf["Y"].values

    z_pred, _ = ok_model.execute(
        "points", grid_x, grid_y,
        backend="loop",
        n_closest_points=10
    )

结果可视化

本节展示土壤属性空间分布图，直观展示插值结果。

有机质空间分布

pH值空间分布

速效钾空间分布

结果解读与结论

插值结果分析

-从空间分布图可以看出，曾都区土壤有机质含量整体呈现南高北低、西南高东北低的空间分异格局：研究区西南部（图中左下角）有机质含量最高，普遍达到 30-35 g/kg；东北部（图中右上角）含量最低，多在 10-15 g/kg 之间，含量梯度自西南向东北呈平滑过渡。这种分布趋势与区域地形、土地利用类型的空间异质性高度契合，也与原始样点的统计特征一致。

-本次普通克里格插值结果整体平滑连续，无明显突变或异常斑块：即使在样点分布稀疏的南部区域，也未出现不合理的数值跳跃，说明模型成功捕捉了土壤有机质的空间自相关性，实现了无偏最优的空间预测，有效避免了传统反距离加权插值易出现的 “牛眼效应”。

-从数值范围来看，插值结果完整覆盖了研究区土壤有机质的变异区间，且含量梯度过渡自然，未出现极端异常值；中北部密集的采样点有效约束了插值精度，南部区域的插值结果也保持了良好的空间连续性，验证了本次插值方案的可靠性。

方法说明

本次仅使用普通克里格（Ordinary Kriging, OK） 进行土壤属性空间插值，该方法是土壤数字制图领域的经典地统计方法，核心优势在于： 1. 能够充分利用样点的空间自相关性，实现无偏最优预测 2. 能够同时输出预测值和误差方差，量化插值结果的不确定性 3. 模型结构简单，参数易于解释，符合县域土壤属性制图的实际应用需求 4. 可在ordinary_kriging_full.py代码中的属性配置中选择不同土壤属性，获取对应数字制图结果

可复现性说明

本项目全部采用 Python 开源工具库开发，通过environment.yml文件即可完成环境配置，没有特殊的环境依赖，通用性强。
代码中的数据读取路径、字段列名都和本次使用的曾都历史数据.xls、曾都渔网点.shp 完全匹配，拿到代码后可以直接运行，不需要额外修改。
本次插值使用的高斯半方差模型、最近 10 个点搜索等核心参数均为固定设置，整个计算过程不包含随机因素，在不同电脑、不同时间运行，都能得到完全一致的结果。