一种基于多重分形的土壤水分特征曲线模拟方法

技术领域

本发明属于农业技术领域，具体说，涉及一种基于多重分形的土壤水分 特征曲线模拟方法。

背景技术

土壤水分特征曲线是描述土壤含水量与吸力(基质势)之间的关系曲线。 它反映了土壤水能量与土壤水含量的函数关系，是进行土壤水分和溶质运移 模拟的重要物理参数，土壤水分特征曲线的模拟精度直接影响了土壤水分和 溶质运移模拟预测的精度，因此它是表示土壤基本水力特性最重要的指标之 一。

土壤水分特征曲线是一个高度非线性的函数，受到多种环境和土壤本身 因素的制约，在空间上和时间上都具有强烈的变异性，准确测定非常困难。 目前确定土壤水分特征曲线的方法主要有直接测定法和间接模拟法两种，直 接方法确定土壤水分特征曲线的过程费时耗力，且测定结果具有较大不确定 性；间接方法由于较少的考虑土壤结构对水分特征曲线模拟的影响，导致土 壤水分模拟的精度不高。

发明内容

本发明的目的在于基于多重分形理论，提供一种间接确定土壤水分特征 曲线的方法，以提高土壤水分特征曲线的模拟精度。

本发明的实施例提供了一种基于多重分形的土壤水分特征曲线模拟方 法，包括：

步骤1，获取土壤数字图像；

步骤2，对该土壤数字图像进行处理，获取土壤孔隙分布特征数据，并 计算土壤孔隙分布的多重分形参数；

步骤3，根据土壤孔隙直径与土壤毛管吸力Young-Laplace方程确定土 壤水分特征曲线中压力水头与土壤孔隙半径之间的关系；

步骤4，建立压力水头与土壤孔隙半径多重分形参数的奇异性指数函数， 并对该函数进行积分变换，获得土壤水分特征曲线模型。

进一步，步骤3具体包括：

假设土壤的孔隙分布分数P_i包含在孔隙当中，并用孔隙半径r_i进行估计：

$P_i=c_1{r}_i^{a_i}---\left(1\right)$

式中c₁为经验常数，α_i为奇异性指数；

孔隙半径r_i、奇异性指数为α_i的孔隙数目N_i之间满足下式：

$N_i=c_2{r}_i^{-f\left({\alpha }_i\right)}---\left(2\right)$

式中c₂为经验常数，f(α_i)为具有奇异性指数α_i的孔隙分形维数，设定多 孔介质的所有孔隙中经验常数c₁和c₂相等，P_i的值只随α_i值变化，则包含在 所有孔隙当中的部分孔隙可以表示为：

$P=N_i{P}_i=c_1{c}_2{r}_i^{{\alpha }_i-f\left({\alpha }_i\right)}---\left(3\right)$

对(3)式积分有：

$P\left(r_i\right)=c_1{c}_2{\int }_{{\alpha }_{\min }}^{{\alpha }_{max}}{r}_i^{\alpha -f\left(\alpha \right)}d\alpha ---\left(4\right)$

式中α_max和α_min分别为奇异性指数α的最大值和最小值；

确定α和f(α)之间的关系，用下面的公式来描述他们之间的关系：

$f\left(\alpha \right)=1-\frac{2}{ln2}\left(\frac{\alpha -{\alpha }_0}{{\alpha }_{max}-{\alpha }_{min}}\right)---\left(5\right)$

式中α₀＝(α_max+α_min)/2；

据土壤毛管吸力Young-Laplace方程获取土壤孔隙半径与压力水头之 间的关系，土壤毛管吸力Young-Laplace方程为：

$h\left(r\right)=\frac{2\sigma cos\omega }{r}---\left(6\right)$

式中：h为压力水头；r为土壤孔隙半径；σ为表面张力；ω为接触 角，在25℃时，水的表面张力σ取0.075N/m，接触角ω取0°，则：

$h=\frac{0.15}{r}---\left(7\right)$

式中h是压力水头。

进一步，步骤4具体包括：

将式(5)、(7)代入式(4):得：

$\begin{array}{c}P\left(h_i\right)=c_1{c}_2\sqrt{\frac{\pi \ln 2}{8\ln \left(h_i\right)}}({\alpha }_{\max }-{\alpha }_{\min }){h_i}^{1/8\left(\ln 2\right){({\alpha }_{\max }-{\alpha }_{\min })}^2-1/2({\alpha }_{\max }+{\alpha }_{\min })+1}\times \\ [erf\left(\sqrt{\frac{\ln \left(h_i\right)}{8\ln 2}}\right(\ln 2{\alpha }_{\min }-\ln 2{\alpha }_{\max }+2\left)\right)-erf\left(\sqrt{\frac{\ln \left(h_i\right)}{8\ln 2}}\right(\ln 2{\alpha }_{\min }-\ln 2{\alpha }_{\max }-2\left)\right)]\end{array}---\left(8\right)$

通过对式(8)的数值积分得到相应于某压力水头h的土壤含水量θ(h)， 即：

$\theta \left(h\right)=c_1{c}_2{\int }_{h_{min}}^h\sqrt{\frac{\pi ln2}{8ln\left(h_i\right)}}({\alpha }_{max}-{\alpha }_{min}){h_i}^{1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2-1/2({\alpha }_{max}+{\alpha }_{min})+1}(a_1-a_2)dh---\left(9\right)$

式中：h_min是相应于最小孔隙的压力水头，a₁和a₂表示为：

$\begin{array}{c}{a}_1=erf\left[\sqrt{\frac{\ln h}{8\ln 2}}(\ln 2{\alpha }_{\min }-ln2{\alpha }_{\max }+2)\right]\\ a_2=erf\left[\sqrt{\frac{\ln h}{8\ln 2}}(\ln 2{\alpha }_{\min }-ln2{\alpha }_{\max }-2)\right]\end{array}---\left(10\right)$

假定土壤孔隙半径r达到最大孔隙半径r_max时的土壤含水量θ(r)等于饱和 含水量θ_s，相应于最大孔隙半径r_max的压力水头h_max为进气口压力，饱和含水 量用下式估计：

${\theta }_s=c_1{c}_2{\int }_{h_{min}}^{h_{max}}\sqrt{\frac{\pi ln2}{8ln\left(h\right)}}({\alpha }_{max}-{\alpha }_{min})h^{1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2-1/2({\alpha }_{max}+{\alpha }_{min})+1}(a_1-a_2)dh---\left(11\right)$

由式(9)和(11)可以得到由参数θ_s，α_max，α_min，h_min和h_max表示的土 壤水分特征曲线模型：

$\frac{\theta \left(h\right)}{{\theta }_s}=\frac{{\int }_{r_{min}}^r\sqrt{\frac{-1}{ln\left(r\right)}}{r}^{-1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2+1/2({\alpha }_{max}+{\alpha }_{min})-1}(a_1-a_2)dr}{{\int }_{r_{min}}^{r_{\max }}\sqrt{\frac{-1}{ln\left(r\right)}}{r}^{-1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2+1/2({\alpha }_{max}+{\alpha }_{min})-1}(a_1-a_2)dr}---\left(12\right)$

式(12)即为基于多重分形的土壤水分特征曲线模型与现有技术相比本 发明的有益效果是：基于多重分形方法，充分考虑考虑土壤结构对水分特征 曲线模拟的影响，提高了土壤水分特征曲线的模拟精度。

附图说明

图1是本发明一种基于多重分形的土壤水分特征曲线模拟方法的流程 图。

具体实施方式

下面结合附图所示的各实施方式对本发明进行详细说明，但应当说明的 是，这些实施方式并非对本发明的限制，本领域普通技术人员根据这些实施 方式所作的功能、方法、或者结构上的等效变换或替代，均属于本发明的保 护范围之内。

参图1所示，图1是本发明一种基于多重分形的土壤水分特征曲线模拟 方法的流程图。

本实施例提供了一种基于多重分形的土壤水分特征曲线模拟方法，包 括：

步骤S1，获取土壤数字图像；

步骤S2，对土壤数字图像进行处理，获取土壤孔隙分布特征数据，并计 算土壤孔隙分布的多重分形参数；

步骤S3，根据土壤孔隙直径与土壤毛管吸力Young-Laplace方程确定土 壤水分特征曲线中压力水头与土壤孔隙半径之间的关系；

步骤S4，建立压力水头与土壤孔隙半径多重分形参数的奇异性指数函数 (将土壤水分特征曲线中的压力水头表示为土壤孔隙多重分形参数的函 数)，并对该函数进行积分变换，获得(基于土壤孔隙分布多重分形奇异性 指数的)土壤水分特征曲线模型。

本实施例提供的土壤水分特征曲线模拟方法，基于多重分形方法，充分 考虑考虑土壤结构对水分特征曲线模拟的影响，提高了土壤水分特征曲线的 模拟精度。

下面详述本发明的具体实施过程：

1)选定研究区域，田间采集土壤样本，通过CT图像扫描获得土壤数字 图像；采用计算机数字图像处理技术和ImageJ图像处理软件获得土壤孔隙分 布状况，并计算土壤孔隙分布的多重分形参数。

2)假设土壤的孔隙分布分数P_i包含在孔隙当中可以用孔隙半径r_i来估 计：

$P_i=c_1{r}_i^{a_i}---\left(1\right)$

式中c₁为经验常数，α_i为奇异性指数，孔隙半径r_i、奇异性指数为α_i的孔 隙数目N_i之间满足下式：

$N_i=c_2{r}_i^{-f\left({\alpha }_i\right)}---\left(2\right)$

式中c₂为经验常数，f(α_i)为具有奇异性指数α_i的孔隙分形维数，我们认 为在所研究的多孔介质的所有孔隙中经验参数c₁和c₂相等，P_i的值只随α_i值 变化，因此包含在所有孔隙当中的部分孔隙可以表示为：

$P=N_i{P}_i=c_1{c}_2{r}_i^{{\alpha }_i-f\left({\alpha }_i\right)}---\left(3\right)$

对(3)式积分有：

$P\left(r_i\right)=c_1{c}_2{\int }_{{\alpha }_{\min }}^{{\alpha }_{max}}{r}_i^{\alpha -f\left(\alpha \right)}d\alpha ---\left(4\right)$

式中α_max和α_min分别为奇异性指数α的最大值和最小值。

3)确定α和f(α)之间的关系，用下面的公式来描述他们之间的关系：

$f\left(\alpha \right)=1-\frac{2}{ln2}\left(\frac{\alpha -{\alpha }_0}{{\alpha }_{max}-{\alpha }_{min}}\right)---\left(5\right)$

式中α₀＝(α_max+α_min)/2。

4)孔隙半径与压力水头之间的关系可以用毛管上升理论方程表示，根 据土壤毛管吸力Young-Laplace方程：

$h\left(r\right)=\frac{2\sigma cos\omega }{r}---\left(6\right)$

式中：h为土壤负压(基质势)；r为孔隙尺寸；σ为表面张力；ω为接 触角。在25℃时，水的表面张力σ约为0.075N/m，一般取接触角ω为0 °。则：

$h=\frac{0.15}{r}---\left(7\right)$

式中h是压力水头。

将式(5)、(7)代入式(4)可得：

$\begin{array}{c}P\left(h_i\right)=c_1{c}_2\sqrt{\frac{\pi \ln 2}{8\ln \left(h_i\right)}}({\alpha }_{\max }-{\alpha }_{\min }){h_i}^{1/8\left(\ln 2\right){({\alpha }_{\max }-{\alpha }_{\min })}^2-1/2({\alpha }_{\max }+{\alpha }_{\min })+1}\times \\ [erf\left(\sqrt{\frac{\ln \left(h_i\right)}{8\ln 2}}\right(\ln 2{\alpha }_{\min }-\ln 2{\alpha }_{\max }+2\left)\right)-erf\left(\sqrt{\frac{\ln \left(h_i\right)}{8\ln 2}}\right(\ln 2{\alpha }_{\min }-\ln 2{\alpha }_{\max }-2\left)\right)]\end{array}---\left(8\right)$

5)相应于某压力水头h的土壤含水量θ(h)可以通过对式(8)的数值积 分得到，即：

$\theta \left(h\right)=c_1{c}_2{\int }_{h_{min}}^h\sqrt{\frac{\pi ln2}{8ln\left(h_i\right)}}({\alpha }_{max}-{\alpha }_{min}){h_i}^{1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2-1/2({\alpha }_{max}+{\alpha }_{min})+1}(a_1-a_2)dh---\left(9\right)$

式中：h_min是相应于最小孔隙的压力水头，a₁和a₂可表示为：

$\begin{array}{c}{a}_1=erf\left[\sqrt{\frac{\ln h}{8\ln 2}}(\ln 2{\alpha }_{\min }-ln2{\alpha }_{\max }+2)\right]\\ a_2=erf\left[\sqrt{\frac{\ln h}{8\ln 2}}(\ln 2{\alpha }_{\min }-ln2{\alpha }_{\max }-2)\right]\end{array}---\left(10\right)$

6)假定土壤孔隙半径r达到最大孔隙半径r_max时的土壤含水量θ(r)等于饱 和含水量θ_s，相应于最大孔隙半径r_max的压力水头h_max为进气口压力，饱和含 水量可以用下式估计：

${\theta }_s=c_1{c}_2{\int }_{h_{min}}^{h_{max}}\sqrt{\frac{\pi ln2}{8ln\left(h\right)}}({\alpha }_{max}-{\alpha }_{min})h^{1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2-1/2({\alpha }_{max}+{\alpha }_{min})+1}(a_1-a_2)dh---\left(11\right)$

由式(9)和(11)可以得到由参数θ_s，α_max，α_min，h_min和h_max表示的土 壤水分特征曲线模型：

$\frac{\theta \left(h\right)}{{\theta }_s}=\frac{{\int }_{r_{min}}^r\sqrt{\frac{-1}{ln\left(r\right)}}{r}^{-1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2+1/2({\alpha }_{max}+{\alpha }_{min})-1}(a_1-a_2)dr}{{\int }_{r_{min}}^{r_{\max }}\sqrt{\frac{-1}{ln\left(r\right)}}{r}^{-1/8\left(ln2\right){({\alpha }_{max}-{\alpha }_{min})}^2+1/2({\alpha }_{max}+{\alpha }_{min})-1}(a_1-a_2)dr}---\left(12\right)$

a₁和a₂可表示为：

$a_1=erf\left[\sqrt{\frac{-ln\left(r\right)}{8ln2}}(ln2{\alpha }_{min}-ln2{\alpha }_{max}-2)\right]$

$a_2=erf\left[\sqrt{\frac{-ln\left(r\right)}{8ln2}}(ln2{\alpha }_{min}-ln2{\alpha }_{max}+2)\right]$

式(12)即为基于多重分形的土壤水分特征曲线模型。

本实施例提供的土壤水分特征曲线模拟方法，采用多重分形方法定量描 述土壤结构的非均质特征，并应用于土壤水分特征曲线的模拟中，提高了土 壤水分特征曲线的模拟精度。

上文所列出的一系列的详细说明仅仅是针对本发明的可行性实施方式 的具体说明，它们并非用以限制本发明的保护范围，凡未脱离本发明技艺精 神所作的等效实施方式或变更均应包含在本发明的保护范围之内。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细 节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形 式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的， 而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因 此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发 明内。