基于van Genuchten模型绘制土壤水分特征曲线的MATLAB问题

你当前的核心问题是误用了幂函数模型（a*(x^b)）进行拟合，而非van Genuchten模型本身，导致拟合结果完全不符合土壤水分特征曲线的规律。以下是问题分析和修正方案：

原代码及拟合结果 #

% Load data from Excel file
EV = readtable('Retention.xlsx', 'Sheet', 1);
CB1 = EV(1:38, :);
theta = CB1.SWC_A / 100;  % to change the data in decimale point
h = CB1.psi;
% Parameters obtained from previous experiments
theta_r = 0.070; % Residual water content
theta_s = 0.434; % Saturated water content
alpha = -3.402; % air-entry value
n = 0.3273; % Pore size distribution parameter
m = 3.4195; % Empirical parameter

%Calculate theta_c using the model equation ((van Genuchten model)
theta_c = theta_r + ((theta_s - theta_r) ./ (1 + alpha * abs(h)).^n).^m;

% Fitting the model to the observed theta values (not theta_c)
fitType = fittype('a * (x^b)');
fitOptions = fitoptions(fitType);
fitOptions.StartPoint = [1 1];
fittedModel = fit(h, theta, fitType, fitOptions);

% Calculate fitted theta values using the fitted model
fittedTheta = fittedModel(h);

%Plotting
figure;
hold on;
scatter(h, theta, 'filled', 'DisplayName', 'Data');
plot(h, fittedTheta, 'r-', 'LineWidth', 1); % Fitted curve
xlabel('Pressure Head (h)');
ylabel('Volumetric Moisture Content (\theta)');
legend('Location', 'SouthEast');
box on;%%%

拟合结果图：

修正方案 #

核心问题说明 #

van Genuchten模型的标准形式为：
$$\theta(h) = \theta_r + \frac{\theta_s - \theta_r}{\left[1 + (\alpha |h|)^n\right]m}$$
通常可约定 $m = 1 - 1/n$（也可作为独立参数拟合），必须直接用该模型对观测数据进行拟合，而非无关的幂函数模型。

修正后的完整代码 #

% 加载Excel数据
EV = readtable('Retention.xlsx', 'Sheet', 1);
CB1 = EV(1:38, :);
theta = CB1.SWC_A / 100;  % 转换为小数格式
h = CB1.psi;

% 定义van Genuchten拟合模型
fitType = fittype(...
    'theta_r + (theta_s - theta_r) ./ (1 + (alpha*abs(h)).^n).^m', ...
    'Independent', 'h', ...
    'Dependent', 'theta');

% 设置拟合参数初始值（基于你已有实验数据）
fitOptions = fitoptions(fitType);
fitOptions.StartPoint = [0.07, 0.434, 3.402, 0.3273, 3.4195]; % theta_r, theta_s, alpha, n, m
% 设置参数物理约束：alpha>0，n>1，含水率在0-1之间
fitOptions.Lower = [0, 0, 1e-3, 1.01, 0];
fitOptions.Upper = [0.2, 0.6, 10, 10, 10];

% 执行拟合
fittedModel = fit(h, theta, fitType, fitOptions);

% 计算拟合值
fittedTheta = fittedModel(h);

% 绘图展示
figure;
hold on;
scatter(h, theta, 'filled', 'DisplayName', '观测数据');
plot(h, fittedTheta, 'r-', 'LineWidth', 1.5, 'DisplayName', 'van Genuchten拟合曲线');
xlabel('压头 (h)');
ylabel('体积含水率 (\theta)');
legend('Location', 'SouthEast');
box on;
grid on;

% 输出拟合得到的参数
disp('拟合得到的van Genuchten参数：');
disp(fittedModel);

关键注意点 #

1. 参数符号：原代码中alpha设为负数，这不符合模型定义（$\alpha$为正的进气值参数），修正时已调整为正的初始值并添加约束。
2. 参数约束：根据模型物理意义设置上下限，避免拟合出无意义的参数值，提升拟合稳定性。
3. 初始值：使用已有实验参数作为拟合初始值，能显著提高拟合的收敛速度和准确性。

内容的提问来源于stack exchange，提问作者Asegidew