【SWAT水文模型】土壤质地转化

1 土壤质地转化

土壤粒径分布 是土壤物理性质的重要指标，主要表示土壤固相中不同粒径颗粒的比例。土壤粒径分布数据通常以累积百分含量曲线表示，并用于估算其他土壤参数，如饱和导水率、有效持水量等。

在不同标准中，土壤粒径分布的分类有所差异：

• 国际制：细分为粘粒、粉砂、细砂、粗砂和石砾。
• 美国制（USDA标准）：简化为粘粒 (CLAY)、粉砂 (SILT)、砂粒 (SAND) 和石砾 (ROCK)。

从中国土壤数据库下载的数据，是我国第二次土壤普查采用的国际制，而SWAT模型采用的土壤粒径级配标准是USDA简化的美制标准，因此，存在一个国际制向美国制转换的问题。

国际制与美国制区别如下表所示：

土壤粒径分布是指土壤固相中不同粗细级别的土粒所占的比例，常用某一粒径及其对应的累积百分含量曲线来表示。

土壤质地转换方法有多种，考虑到模型的通用性，参数形式的土壤粒径分布模型更便于标准程序的编制以及不同来源粒径分析资料的对比和统一。

2 转换模型

2.1 双参数修正的经验逻辑生长模型

2.1.1 原理

采用双参数修正的经验逻辑生长模型进行土壤粒径分布的标准化转换，其模型定义如下：此模型将每类土壤质地的粒径完全累积分布看作1

式中，P(d)是土壤颗粒直径小于d 的累积分布；d0 是模型运行的最小土壤颗粒直径，u、c 是模型参数。

参数作用：

• 参数 $u$：控制分布曲线的陡峭程度；
• 参数 $c$：调整分布的非对称性。

注意事项
1、参数调整：根据不同地区的实际土壤特性，参数 $u$ 和 $c$ 需要调整以优化拟合效果。

2、输入数据格式：确保输入的土壤粒径分布数据符合国际制分类。

3、模型假设：假设土壤质地的粒径分布遵循经验逻辑生长模型，可能存在一定误差。

2.1.2 代码实现（MATLAB）

以下为基于经验逻辑生长模型的土壤质地转换代码实现：

% 土壤粒径分布转换模型 - MATLAB实现
% 根据土壤粒径分布模型实现国际制向美国制转换

%% 清空环境变量
clear; 
clc;

% 模型参数定义
u = 2;      % 调整陡峭程度的参数
c = 1.5;    % 调整非对称性的参数
d0 = 0.002; % 最小土壤颗粒直径（单位：mm）

% 粒径范围（国际制分类）
d_international = [0.002, 0.02, 0.2, 2]; % 国际制粒径边界

% 累积分布函数 P(d) 的前置计算
% P(d0) 是初始累积分布值，假设为 0.1（可以根据实际情况调整）
P_d0 = 0.1;

% 定义累积分布函数 P(d)
P = @(d) 1 ./ (1 + (1 / P_d0 - 1) .* exp(-u .* ((d - d0) ./ d0).^c));

% 定义粒径范围内的均匀分布（用于绘制分布曲线）
d_range = linspace(0, 2, 1000); % 从 0 mm 到 2 mm 的粒径范围

% 计算每个粒径的累积百分含量
P_values = P(d_range);

% 绘制土壤粒径累计分布曲线
figure;
plot(d_range, P_values, 'LineWidth', 2);
xlabel('粒径 d (mm)');
ylabel('累积百分含量 P(d)');
title('土壤粒径累积分布曲线');
grid on;

%% 将国际制转换为美国制
% 国际制类别：黏粒、粉砂、细砂、粗砂、石砾
% 美国制类别：CLAY, SILT, SAND, ROCK
clay_percent = P(0.002); % 粘粒（<0.002mm）
silt_percent = P(0.05) - clay_percent; % 粉砂（0.002-0.05mm）
sand_percent = P(2) - (silt_percent + clay_percent); % 砂粒（0.05-2mm）
rock_percent = 1 - P(2); % 石砾（>2mm）

% 输出转换结果
fprintf('美国制土壤粒径百分比分布结果：\n');
fprintf('黏粒 (CLAY): %.2f%%\n', clay_percent * 100);
fprintf('粉砂 (SILT): %.2f%%\n', silt_percent * 100);
fprintf('砂粒 (SAND): %.2f%%\n', sand_percent * 100);
fprintf('石砾 (ROCK): %.2f%%\n', rock_percent * 100);

% 可视化国际制与美国制转换对比
categories = {'CLAY', 'SILT', 'SAND', 'ROCK'};
percentages = [clay_percent, silt_percent, sand_percent, rock_percent] * 100;

figure;
bar(categorical(categories), percentages, 'FaceColor', [0.2, 0.6, 0.8]);
xlabel('土壤类别');
ylabel('累积百分含量 (%)');
title('国际制向美国制土壤质地转换结果');
grid on;

参考