用于预测根区土壤水分分布的LRBF无网格方法（Matlab代码实现）-优快云博客

本文链接：https://blog.youkuaiyun.com/F_Matlab/article/details/148931060

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

⛳️赠与读者

💥1 概述

一、LRBF无网格方法的基本原理与定义

二、根区土壤水分分布的关键影响因素

1. 根系特征（直接主导水分空间分布）

2. 土壤特性（通过改变结构间接调控水分）

3. 外部环境（动态边界条件）

三、LRBF方法在根区水分预测中的应用案例

案例1：蒸发过程模拟（1D）

案例2：根系吸水动态（1D）

案例3：沟灌与圆灌模拟（2D/3D）

四、LRBF与传统网格方法的对比分析

五、预测精度与参数敏感性的研究进展

1. 精度验证

2. 参数敏感性

3. 优化方向

六、挑战与未来展望

结论

📚2 运行结果

2.1 算例1：蒸发过程中的土壤水分动态

2.2 算例2：变表面通量条件下根系土壤中的非饱和流

2.3 算例3：周期性灌溉沟渠的二维水流

2.4 算例4：圆形水源的3D灌溉

🎉3 参考文献

🌈4 Matlab代码、数据、文章下载

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

在本文中，我们使用空间中的局部径向基函数（LRBF）无网格技术和时间离散化的向后欧拉方案来求解该系统。LRBF方法是一种精确且计算效率高的方法，不需要网格生成，在解决高维问题时具有灵活性。此外，这种方法产生了一个稀疏矩阵系统，避免了病态问题。该数值模型实现了对一维、二维和三维土壤的渗透和植物根系吸水。使用非平凡解析解和可用实验数据进行数值实验，以验证所提出的数值技术的性能。

在本文中，我们首先提出了土壤中非饱和流动和植物根系吸水的耦合数值模型。在开发的数值模型中使用了理查兹方程和不同的公式来描述根区的渗透，并研究植物根系对土壤水分分布的影响。使用基尔霍夫变换的理查兹方程，并考虑加德纳模型用于毛细管压力。在我们的方法中，我们采用基于局部径向基函数（LRBF）的无网格方法来求解得到的方程组。LRBF方法是一种精确且计算效率高的方法，不需要网格生成，在处理具有复杂几何形状的高维问题时具有灵活性。此外，这种方法产生了一个稀疏矩阵系统，避免了病态问题。我们实现了渗透和植物根系吸水的耦合数值模型，用于一维、二维和三维土壤。使用非平凡解析解和可用实验数据进行数值实验，以验证耦合数值模型。数值结果表明了所提出的数值方法在预测根区土壤水分动态方面的性能和能力。

在本文中，我们重点对根区土壤水分分布进行建模引入计算技术，以有效地解决耦合数值模型土壤渗透和植物根系吸水。所提出的模型基于理查兹方程和不同的根系吸水公式来研究影响植物根系吸水对土壤水分分布的影响。加德纳模型用于毛细管压力。在系统中加入根系吸水的隐含汇项在数值分辨率方面增加了复杂性。基尔霍夫变换采用控制方程来降低系统的非线性。我们的数字该方法基于LRBF无网格法求解耦合数值模型。这个该技术提供了显著的计算优势，包括减少计算时间由于不需要网格生成，因此它具有成本低和数值解精确的优点。我们提出的耦合模型通过各种数值试验进行了验证，包括一、二、三、四层土壤在蒸发和根系吸水过程中的土壤水分分布，以及三维案例。我们的方法的可靠性得到了实验数据和非平凡精确解的验证。数值结果表明了该算法的性能
有效预测土壤中非饱和流动态的拟议方法在蒸发和植物根系吸水的作用下。我们的研究提供了一个稳健的数值研究根系吸水对土壤水分动态影响的框架，即对农业和生态系统管理很重要。拟议的数值技术可以进一步扩展，以纳入更复杂的根系吸水模型可以为未来的研究提供有前景的途径。

一、LRBF无网格方法的基本原理与定义

LRBF（Local Radial Basis Function）无网格法是一种基于径向基函数插值的数值模拟技术。其核心思想是在求解域内布置离散节点，每个节点周围定义局部支持域，利用支持域内节点的函数值构建局部逼近函数求解偏微分方程（如Richards方程）。与传统网格法（如有限元、有限差分）相比，其关键创新点包括：

节点离散化：无需预划分网格，节点可规则或非规则分布，对节点分布的敏感性较低。
局部支持域：每个节点仅依赖邻近节点信息，避免全局矩阵求解，提升计算效率。
Richards方程求解：通过基尔霍夫变换（Kirchhoff Transform）降低方程非线性，结合后向欧拉法（BDF）进行时间离散化和Picard迭代线性化系统。

二、根区土壤水分分布的关键影响因素

LRBF方法需耦合以下影响因子以提升预测精度：

1. 根系特征（直接主导水分空间分布）

根系生物量（RB）与根长密度（RLD） ：与染色面积比（DAR）、土壤饱和含水率（SMC）呈显著正相关（R²≈0.45）。毛竹林因RB和RLD最高（休耕地的30倍），其60cm土层水分分布均匀性显著优于休耕地。
根系形态：比根长对水分湍动强度影响最大（灰色关联度0.984），粗根（>1mm）促进水分下渗，细根（<0.5mm）则限制水分至浅层（0–20cm）。

2. 土壤特性（通过改变结构间接调控水分）

团聚体稳定性：根系分泌物促进土壤颗粒聚集，根系生物量（RB）对团聚体稳定性（MWD）的影响（R²=0.83）显著大于有机质（R²=0.38）。
非均质性：黏重土壤或高地下水位限制根系分布；肥沃土壤中水肥气热协调促进根系富集。

3. 外部环境（动态边界条件）

气候脉冲：降水时间、强度与频率直接影响根土界面水分再分配。
土壤水势梯度：土壤含水量降低时，根冠比增大以提升吸水效率，但过度庞大的根系消耗光合产物，反降低产量。

三、LRBF方法在根区水分预测中的应用案例

案例1：蒸发过程模拟（1D）

场景：模拟K-7砂土在3种温湿度组合下的蒸发实验。
方法：LRBF求解Richards方程，对比实测数据与解析解。
结果：成功捕捉蒸发两阶段（恒速期90–180小时，降速期），含水量预测RMSE<0.05。

案例2：根系吸水动态（1D）

场景：100cm深根土，地表通量变化（稳态/瞬态）耦合两种根系吸水模型（指数型/阶梯型）。
结果：
- 根系吸水显著缩短稳态达时间（无根系需50小时，有根系仅30小时）。
- 瞬态通量下，LRBF精确量化根区-底土界面水通量，RMSE趋近于0。

案例3：沟灌与圆灌模拟（2D/3D）

沟灌（2D） ：模拟周期性格栅沟入渗，LRBF处理非均质土壤通量边界，解析水分横向扩散。
圆灌（3D） ：圆柱域内圆形灌溉源模拟，LRBF耦合贝塞尔函数解析解，验证吸水参数对湿润区半径的影响。

四、LRBF与传统网格方法的对比分析

维度	LRBF无网格法	传统网格法（有限元/差分）
前处理	无需网格生成，节点布置灵活	复杂几何网格划分耗时，易畸变
非均质性处理	节点级参数定义，直接区分土壤类型	需细分网格区分材质，增加计算负担
根系吸水模拟	在根系节点直接施加吸水项，精度高	根系影响均摊至网格单元，局部失真
复杂边界	自然处理动态通量/根系边界	网格重构代价高，易数值振荡
高维问题	易于扩展至3D模拟	网格数量指数增长，效率骤降
优势总结：LRBF显著简化前处理，提升非均质环境和离散根系模拟精度。

五、预测精度与参数敏感性的研究进展

1. 精度验证

实验数据对比：在蒸发、入渗、根系吸水案例中，LRBF预测值与实验数据/解析解的均方根误差（RMSE）均低于0.1，证实其可靠性。
高维扩展能力：3D灌溉模拟中，LRBF生成的稀疏矩阵避免病态问题，计算效率优于全局RBF。

2. 参数敏感性

形状参数（c） ：径向基函数宽度影响局部逼近精度。研究表明自适应c选择策略可平衡精度与稳定性。
支持域节点数（ns） ：过少导致欠拟合，过多增加计算量。推荐ns=5–20，结合KD树加速邻近节点搜索。
停止容差（stopTol） ：Picard迭代容差设为1e-5可确保收敛且避免过度计算。

3. 优化方向

机器学习辅助：RBF神经网络可修正模型参数（如导管动力学预测误差降低75%）。
不确定性量化：需整合土壤参数、根系分布的先验分布，通过蒙特卡洛提升预测鲁棒性。

六、挑战与未来展望

算法优化：开发形状参数自适应选择策略，结合SPH/MLS等多方法优势。
生理过程耦合：集成根系水力导度模型与植物蒸腾反馈机制，提升吸水项动态精度。
数据-模型融合：同化遥感土壤水分数据（如SMAP）与地面传感器网络，校准空间异质性。
生态水文扩展：链接根区水分预测与边坡稳定性模型，量化森林根系的水土保持效应。

结论

LRBF无网格法通过节点离散化和局部支持域构建，克服了传统网格法在根区复杂环境（非均质土壤、离散根系、动态边界）中的局限性，显著提升Richards方程的求解效率和精度。未来需重点突破参数自适应优化、多过程耦合及不确定性量化，结合遥感与地面观测数据，推动根区水分模拟向高分辨率、动态预报方向发展。

📚2 运行结果

2.1 算例1：蒸发过程中的土壤水分动态

2.2 算例2：变表面通量条件下根系土壤中的非饱和流

2.3 算例3：周期性灌溉沟渠的二维水流

2.4 算例4：圆形水源的3D灌溉

部分代码：

% =========================================================================
% Computational Domain Setup
% =========================================================================
xinf = 0; xsup = l; yinf = 0; ysup = 8; % Domain limits [cm]
ny = 1500; % Number of grid points in the y direction
nx = floor(xsup * ny / ysup); % Number of grid points in the x direction
dx = (xsup - xinf) / (nx - 1); % Grid spacing in the x direction [cm]
dy = (ysup - yinf) / (ny - 1); % Grid spacing in the y direction [cm]
n = nx * ny; % Total number of grid points
ns = 5; % Number of nearest neighbors for RBF
c = 0.3; % Shape parameter for RBF

% Generate the computational grid
x = xinf:dx:xsup;
z = yinf:dy:ysup;
[xx, zz] = meshgrid(z, x); % Create a grid of points
X = [zz(:) xx(:) 5*ones(n, 1)]; % Initialize grid points with a third column of 5

% =========================================================================
% Define Time Variables
% =========================================================================
dt = 0.01; % Time step [h]
tmin = 0; % Initial time [h]
T = ts; % Final time [h]
t = (tmin:dt:T); % Time vector [h]
nt = length(t); % Number of time steps

% =========================================================================
% Boundary Condition Setup
% =========================================================================
% Define boundary conditions based on the position in the grid
ind = (X(:, 2) == ysup); % Top boundary
X(ind, 3) = 2;
ind = (X(:, 2) == yinf & 0 <= X(:, 1) & X(:, 1) < x0); % Bottom left boundary
X(ind, 3) = 3;
ind = (X(:, 2) == yinf & x0 < X(:, 1) & X(:, 1) <= l); % Bottom right boundary
X(ind, 3) = 4;
ind = (X(:, 1) == xinf | X(:, 1) == xsup); % Side boundaries
X(ind, 3) = 1;

% Extract internal and boundary nodes
ind = (X(:, 3) == 5); % Internal nodes
X1 = X(ind, 1:2);
g = length(X1);

ind = (X(:, 3) == 1); % Side boundary nodes
X1 = [X1; X(ind, 1:2)];
g1 = length(X1) - g;

ind = (X(:, 3) == 2); % Top boundary nodes
X1 = [X1; X(ind, 1:2)];
g2 = length(X1) - (g + g1);

ind = (X(:, 3) == 3); % Bottom left boundary nodes
X1 = [X1; X(ind, 1:2)];
g3 = length(X1) - (g + g1 + g2);

ind = (X(:, 3) == 4); % Bottom right boundary nodes
X1 = [X1; X(ind, 1:2)];
g4 = length(X1) - (g + g1 + g2 + g3);

% =========================================================================
% KD-Tree Initialization and Memory Allocation
% =========================================================================
tree = kdtree_build(X1); % Build KD-tree for nearest neighbors search
o = ones(1, ns); % Vector of ones for convenience
A1 = zeros(n, 1); Lin = zeros(n, 1); Col = zeros(n, 1); B = zeros(n, 1);

% =========================================================================
% Initial Conditions
% =========================================================================
theta0 = Exact1(X1(:, 1), X1(:, 2), 0); % Calculate exact initial condition
Phi0 = PHI(theta0); % Initial potential
Thetaa = zeros(n, nt); Thetaa(:, 1) = theta0; % Initialize theta matrix
HPHI = zeros(n, nt); HPHI(:, 1) = Phi0; % Initialize potential matrix

Ph1 = Phi0; % Initial guess for the iterative solver
kt1 = 1;

% =========================================================================
% Iterative Time-Stepping Solution
% =========================================================================
stopTol = 1e-5; % Stopping tolerance for Picard iteration

for kt = 2:nt
kt % Display current time step
Phin = HPHI(:, kt-1); % Previous potential
Phimp1n = Ph1; % Initial guess for current potential

psimp1n = Theta(Phimp1n); % Calculate corresponding theta
stop_flag = 0;

% Picard Linearization
while stop_flag == 0
Dmp1n = D(psimp1n); % Compute D(theta)
Rmp1n = R(psimp1n); % Compute R(theta)

% Compute RBF coefficients for each region
k = 1; kk = 1:ns;

% Internal nodes
for i = 1:g
indic = kdtree_k_nearest_neighbors(tree, X1(i,:), ns);
xc = X1(indic, 1); yc = X1(indic, 2);
rx = xc*o - (xc*o)'; ry = yc*o - (yc*o)';
P = RBF(rx, ry, c);
rx = X1(i, 1) - xc; ry = X1(i, 2) - yc;
H = RBF(rx, ry, c) ./ (dt * Dmp1n(i)) - DivRBF(rx, ry, c) + DZRBF(rx, ry, c);
AA = H' / P;
A1(k:k+ns-1) = AA(kk); Lin(k:k+ns-1) = i*o'; Col(k:k+ns-1) = indic;
k = k + ns;
end

% Side boundary nodes
for i = g+1:g+g1
indic = kdtree_k_nearest_neighbors(tree, X1(i,:), ns);
xc = X1(indic, 1); yc = X1(indic, 2);
rx = xc*o - (xc*o)'; ry = yc*o - (yc*o)';
P = RBF(rx, ry, c);
rx = X1(i, 1) - xc; ry = X1(i, 2) - yc;
H = DXRBF(rx, ry, c);
AA = H' / P;
A1(k:k+ns-1) = AA(kk); Lin(k:k+ns-1) = i*o'; Col(k:k+ns-1) = indic;
k = k + ns;
end

% Top boundary nodes
for i = g+g1+1:n
indic = kdtree_k_nearest_neighbors(tree, X1(i,:), ns);
xc = X1(indic, 1); yc = X1(indic, 2);
rx = xc*o - (xc*o)'; ry = yc*o - (yc*o)';
P = RBF(rx, ry, c);
rx = X1(i, 1) - xc; ry = X1(i, 2) - yc;
H = RBF(rx, ry, c) - DZRBF(rx, ry, c);
AA = H' / P;
A1(k:k+ns-1) = AA(kk); Lin(k:k+ns-1) = i*o'; Col(k:k+ns-1) = indic;
k = k + ns;
end

% Construct the right-hand side vector B
B(1:g) = Phin(1:g) ./ (dt * Dmp1n(1:g)) + Rmp1n(1:g);
B(g+1:g+g1) = 0;
B(g+g1+1:g+g1+g2) = exp(A * t(kt)) * F0 * x0 / l * exp(-X1(g+g1+1:g+g1+g2, 2) / 2 * (sqrt(1 - 4 * k1) - 1));
B(g+g1+g2+1:g+g1+g2+g3) = exp(A * t(kt)) * F0;
B(g+g1+g2+g3+1:n) = 0;

% Solve the linear system
SA = sparse(Lin, Col, A1);
Phinp1mp1 = SA \ B;

% Check convergence
deltam = Phinp1mp1 - Phimp1n;
max(abs(deltam(1:n))) % Display the maximum change for convergence check
if max(abs(deltam(1:n))) < stopTol
stop_flag = 1;
Phinp1mp1(1:n) = Phimp1n(1:n) + deltam(1:n);
else
Phinp1mp1(1:n) = Phimp1n(1:n) + deltam(1:n);
Phimp1n = Phinp1mp1;
psimp1n = Theta(Phimp1n);
end
end

Ph1 = Phinp1mp1;
kt1 = kt;
HPHI(:, kt1) = Phinp1mp1;
Thetaa(:, kt1) = Theta(Phinp1mp1);
end

% =========================================================================
% Exact Comparaison and Error Calculation
% =========================================================================
E1 = Exact1(X1(:, 1), X1(:, 2), T); % Compute exact solution at final time
diff = E1 - Thetaa(:, nt); % Difference between exact and computed solutions

% Compute errors
MAE = max(abs(diff)); % Maximum Absolute Error
RMSE = norm(diff) / sqrt(n); % Root Mean Squared Error
RelER = norm(abs(diff), 1) / norm(E1, 1); % Relative Error

% Display errors

RMSE % Display the Root Mean Squared Error

% =========================================================================
% Visualization
% =========================================================================
DT = delaunay(X1(:, 1), X1(:, 2)); % Delaunay triangulation

% Plot the numerical solution
figure(1);
trisurf(DT, X1(:, 1), X1(:, 2), Thetaa(:, nt) ./ thetas);
view(0, -90);
set(gca, 'xaxisLocation', 'top');
shading interp;
colormap jet;

% Plot the exact solution
figure(2);
trisurf(DT, X1(:, 1), X1(:, 2), E1 ./ thetas);
view(0, -90);
set(gca, 'xaxisLocation', 'top');
shading interp;
colormap jet;

figure(3);
trisurf(DT, X1(:, 1), X1(:, 2), Thetaa(:, nt) ./ thetas);
view(0, -90);
set(gca, 'xaxisLocation', 'top');
shading interp;
colormap jet;
ylim([0 1.5]); % Limit the y-axis to the range [0, 1.5]
% Plot the exact solution
figure(4);
trisurf(DT, X1(:, 1), X1(:, 2), E1 ./ thetas);
view(0, -90);
set(gca, 'xaxisLocation', 'top');
shading interp;
colormap jet;
ylim([0 1.5]); % Limit the y-axis to the range [0, 1.5]
% =========================================================================
% Contour Plots for Comparisons
% =========================================================================
kt2 = nt;
C1 = Thetaa(g+1:2:g+2*ny, kt2);
C2 = Thetaa(g+2:2:g+2*ny, kt2);
L1 = Thetaa(g+g1+1:g+g1+g2, kt2);
L2 = Thetaa(g+g1+g2+1:n, kt2-1);
B2 = reshape(Thetaa(1:g, kt2), [nx-2, ny-2]);
Matrix1 = [L2'; B2'; L1'];
MatrixSol = [C1, Matrix1, C2];
Thet = MatrixSol ./ thetas;

% Exact solution matrix
C11 = E1(g+1:2:g+2*ny);
C21 = E1(g+2:2:g+2*ny);
L11 = E1(g+g1+1:g+g1+g2);
L21 = E1(g+g1+g2+1:n);
B21 = reshape(E1(1:g), [nx-2, ny-2]);
Matrix11 = [L21'; B21'; L11'];
MatrixSol = [C11, Matrix11, C21];
Ex = MatrixSol ./ thetas;

num_levels = 5; % Number of contour levels

% Plot contours for z=8

figure(5);
hold on;
contour(x, z, Thet, num_levels, 'ShowText', 'on');
xlabel('x^{*} = x/l_s');
ylabel('z^{*} = z/l_s');
hold on;
contour(x, z, Ex, num_levels);
colormap('parula');
view(0, -90);
set(gca, 'xaxisLocation', 'top');

% Plot contours just for z=1.5

num_levels1 = 25; % Number of contour levels
figure(6);
hold on;
contour(x, z, Thet, num_levels1, 'ShowText', 'on');
xlabel('x^{*} = x/l_s');
ylabel('z^{*} = z/l_s');
hold on;
contour(x, z, Ex, num_levels1);
colormap('parula');