CN114117766B - 一种针对vg模型的参数拟合方法 - Google Patents

Abstract

本发明公开了一种针对VG模型的参数拟合方法。属于土水特征曲线领域，具体操作步骤：分类选取不同的初值对已知土体的吸力和含水率的数据进行拟合，将拟合得到的饱和含水率替换原始初值；将得到的饱和含水率代入初值，将拟合得到的经验参数n和经验参数α进行调整后设为初值，对已知数据进行拟合；将得到的饱和含水率、经验参数n、经验参数α和残余含水率进行调整后设为初值，对已知数据进行拟合；将得到的拟合数值和原始值进行对比，通过异常值检验公式，对原始数据中的异常值进行识别。本发明基于上述算法确定的参数初始值，可使拟合更加高效，得到的参数值更接近真实值。本发明为土水特征曲线的拟合提供参考，使更高效、精确地进行参数优化。

Description

一种针对VG模型的参数拟合方法

技术领域

本发明属于非饱和土力学中的土水特征曲线领域，涉及一种针对VG模型的参数拟合方法，具体的是，涉及一种针对VG模型的拟合初值确定方法。

背景技术

土壤水土特征曲线是反映包气带基质势和土壤含水率之间关系的曲线，以此表示包气带中水分的能量与数量之间的关系；为了计算方便，常将数量关系拟合为经验公式；Van Genuchten模型是由Van Genuchten在1980年提出，是目前采用较多的压力水头与含水率之间的经验关系式之一；VG模型如下所示：

式中，θ为体积含水率；θ_r为残余含水率；θ_s为饱和含水率；h为压力水头；α和n为未知参数。通过对试验获得的压力水头和含水率数据进行拟合，得到Van Genuchten模型中每个参数的值。

现有技术中；采用拟合软件进行参数求解的过程中需要输入估计的初始值，以此初始值为依据进行计算拟合；VG模型中存在饱和含水率θ_s、残余含水率θ_r、经验参数n和经验参数α；在没有任何经验的情况下，想要快速准确地拟合出四个参数，显得十分困难；目前，对于土水特征曲线参数的拟合，大多在于阐述使用不同软件对于土水特征曲线的拟合，而很少有人研究初值对于拟合的影响，以及如何设置初值从而达到快速准确拟合的效果；在拟合过程中，输入四个参数的初始值与参数拟合过程直接相关，当输入的初始值与真实值有较大偏离时，拟合参数根本无法准确反映真实情况，甚至出现拟合失败。

发明内容

发明目的：本发明目的是基于对VG模型深入研究的基础上，提出了VG模型拟合初值设置的一种算法和异常值检验方法。

技术方案：本发明所述的一种针对VG模型的参数拟合方法，具体实现步骤如下：

(1)、根据不同土的分类，选取VG模型中饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r四个参数的初值，对已知土体的吸力水头h和体积含水率θ的数据进行拟合，从而得到饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r的拟合值，

将初值中的饱和含水率θ_S用得到的饱和含水率θ_S拟合值进行替换，保持其他三个参数的初值不变，后再进行拟合；直至第N+1次的饱和含水率拟合值与第N次的饱和含水率拟合值之差的绝对值小于0.0001，且第N+1次的拟合优度与第N次的拟合优度之差的绝对值小于0.0001；

(2)、将步骤(1)中最终进行拟合得到的饱和含水率θ_S、经验参数n和经验参数α的拟合值替换步骤(1)中选取的三个对应参数的初值，保持残余含水率θ_r的初值不变，对已知数据进行拟合，直至第M+1次的经验参数n的拟合值与第M次的经验参数n的拟合值之差的绝对值小于0.002，直至第M+1次的经验参数α的拟合值与第M次的经验参数α的拟合值之差的绝对值小于0.0000001，且第M+1次的拟合优度与第M次的拟合优度之差的绝对值小于0.00001；

(3)、将步骤(2)最终进行拟合得到的饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r的拟合值替换步骤(1)中选取的四个对应参数的初值，对已知数据进行拟合，直至第i+1次的残余含水率θ_r的拟合值与第i次的残余含水率θ_r的拟合值之差的绝对值小于0.0001，且第i+1次的拟合优度与第i次的拟合优度之差的绝对值小于0.000001；

(4)、将步骤(3)最终进行拟合得到的拟合数值和原始值进行对比，通过异常值检验公式，对原始数据中的异常值进行识别。

进一步的，在步骤(1)中，所述不同土的分类共分为四大类，包括普通土体、砂土类、粘土类和混凝土类；

所述普通土体包括粘壤土、壤土、粉砂及粉壤土；

所述砂土类包括壤质砂土和和砂；

所述粘土类包括粉质粘土、粉质粘壤土和砂质粘土；

所述混凝土包括塑性混凝土和混凝土。

进一步的，在步骤(1)中，所述土体的饱和含水率θ_S初值根据低转速下的体积含水率确定，取最低转速下的体积含水率为饱和含水率θ_S的初值、且保留四位有效数字；

其中，所述普通土体的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.01、1.4和0.001；

所述砂土类的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.01、2.2和0.001；

所述粘土类的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.01、1.1和0.01；

所述混凝土的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.001、1.5和0.00001。

进一步的，在步骤(2)中，所述的经验参数n保留四位有效数字设置为下一次拟合初值；所述经验参数α保留三位有效数字设置为下一次拟合初值。

进一步的，在步骤(3)中，所述残余含水率θ_r保留一位有效数字设置为下一次拟合初值。

进一步的，在步骤(1)、(2)及(3)中，所使用的VG模型如下所示：

其中，θ为体积含水率；h为压力水头；θ_r为残余含水率；θ_s为饱和含水率；α和n为经验参数；θ和h为已知参数；θ_r、θ_S、α和n为待定参数。

进一步的，在步骤(1)(2)及(3)中，针对VG模型，使用MATLAB中lsqcurvefit的拟合方法对已知数据进行拟合；具体设置包括初值的设置、待定参数上下限的设置、lsqcurvefit拟合函数的使用，具体的拟合代码如下：

param0＝[0.00001,1.5,0.001,0.3050]；

lb＝[0,1,0,0]；

ub＝[5,5,1,1]；

param＝lsqcurvefit('curvefun',param0,h_data,theta_data,lb,ub,options)

param0为初值设置；中括号内按照前往后的顺序，依次为VG模型中经验参数α、经验参数n、残余含水率θ_r和饱和含水率θ_S；

lb和ub分别为下限设置和上限设置。中括号内按照前往后的顺序，依次为经验参数α、经验参数n、残余含水率θ_r和饱和含水率θ_S；

param为经验参数α、经验参数n、残余含水率θ_r和饱和含水率θ_S的拟合值；

lsqcurvefit为MATLAB中基于最小二乘法的内置拟合方法；

curvefun为待拟合的具体函数，此处函数为VG模型。

进一步的，在步骤(1)(2)及(3)中，所述拟合优度用于度量拟合曲线对于原始数据拟合效果的好坏，其公式如下所示：

其中，R²表示拟合优度，y表示待拟合的数据，表示y的均值，/>表示拟合数据。

进一步的，在步骤(4)中，所述的异常值检验公式如下：

式中，e_i为修正误差贡献率，θ_i为第i次的体积含水率，为经过第一步至第三步拟合的第i次的预测体积含水率，η_i为第i次的修正系数；

所述的修正系数计算公式为式中，h为压力水头，n为拟合最优的经验参数；

当e_i大于0.1时，表示第i次原始数据存在异常，可再次试验进行验证。

有益效果：本发明与现有技术相比，本发明的特点:本发明基于对VG模型深入研究基础上，提出了VG模型拟合初值设置的一种算法和异常值检验方法，首先确定饱和含水率的范围使得后期拟合过程平稳，继而联合确定对于拟合影响最大的经验参数α和经验参数n，最后确定残余含水率，由此确定了不同土体VG模型的参数初始值；上述算法是一种系统且逐步稳健靠近真实值的拟合方案，基于上述算法确定的参数初始值，可以使拟合更加高效，得到的参数值更接近真实值；本发明为土水特征曲线的拟合提供参考，使更高效、精确地进行参数优化。

附图说明

图1是本发明的操作流程图；

图2是本发明实施例中编号1土样的示意图；

图3是本发明实施例中编号2土样的示意图。

具体实施方式

以下结合附图和具体实施例，对本发明做出进一步说明。

如图所述，本发明所述的一种针对VG模型的参数拟合方法，具体实现步骤如下：

第一步：根据不同土的分类，选取VG模型中饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r四个参数的初值，对已知土体的吸力水头h和体积含水率θ的数据进行拟合，从而得到饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r的拟合值，

将初值中的饱和含水率θ_S用得到的饱和含水率θ_S拟合值进行替换，保持其他三个参数的初值不变，后再进行拟合；直至第N+1次的饱和含水率拟合值与第N次的饱和含水率拟合值之差的绝对值小于0.0001，且第N+1次的拟合优度与第N次的拟合优度之差的绝对值小于0.0001；

其中，所述的土体分类，共分为四大类，包括普通土体、砂土类、粘土类和混凝土类；其中，普通土体包括粘壤土、壤土、粉砂、粉壤土等；砂土类包括壤质砂土和和砂；粘土类包括粉质粘土、粉质粘壤土和砂质粘土；混凝土包括塑性混凝土和混凝土；

所述的土体的饱和含水率初值根据低转速下的体积含水率确定，取最低转速下的体积含水率为饱和含水率的初值，且保留四位有效数字；

所述的普通土体的残余含水率、经验参数n和经验参数α的初值分别为0.01、1.4和0.001；

所述的粘土类的残余含水率、经验参数n和经验参数α的初值分别为0.01、1.1和0.01；

所述的砂土类的残余含水率、经验参数n和经验参数α的初值分别为0.01、2.2和0.001；

所述的混凝土的残余含水率、经验参数n和经验参数α的初值分别为0.001、1.5和0.00001。

第二步：将步骤(1)中最终进行拟合得到的饱和含水率θ_S、经验参数n和经验参数α的拟合值替换步骤(1)中选取的三个对应参数的初值，保持残余含水率θ_r的初值不变，对已知数据进行拟合，直至第M+1次的经验参数n的拟合值与第M次的经验参数n的拟合值之差的绝对值小于0.002，直至第M+1次的经验参数α的拟合值与第M次的经验参数α的拟合值之差的绝对值小于0.0000001，且第M+1次的拟合优度与第M次的拟合优度之差的绝对值小于0.00001；

所述的经验参数n和经验参数α在下次拟合前必须同时对原始初值进行替换，保持一致性；

所述的经验参数n保留四位有效数字设置为下一次拟合初值，经验参数α保留三位有效数字设置为下一次拟合初值；

所述的经验参数n和经验参数α拟合至最优，指的是满足此条件下(对于混凝土类的土体：|n_i-n_i-1|≤0.001和|α_i-α_i-1|≤0.00001；其他的土体：|n_i-n_i-1|≤0.01和|α_i-α_i-1|≤0.001)，拟合优度R²在饱和含水率、经验参数n和经验参数α基本不变的条件下到达最大值；

第三步：将步骤(2)最终进行拟合得到的饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r的拟合值替换步骤(1)中选取的四个对应参数的初值，对已知数据进行拟合，直至第i+1次的残余含水率θ_r的拟合值与第i次的残余含水率θ_r的拟合值之差的绝对值小于0.0001，且第i+1次的拟合优度与第i次的拟合优度之差的绝对值小于0.000001；

所述的残余含水率保留两位有效数字设置为下一次拟合初值；

所述的残余含水率拟合达到最优，指的是对于所有土体：|Q_r,i-Q_r,i-1|≤0.01

拟合优度R2在饱和含水率、经验参数n、经验参数α和残余含水率基本不变的条件下到达最大值；

上述的拟合函数为lsqcurvefit函数，拟合代码如下A＝1+(abs(alpha*h)).^n，theta＝theta_r+(theta_s-theta_r)./A.^(1-1/n)；

第四步：将步骤(3)最终进行拟合得到的拟合数值和原始值进行对比，通过异常值检验公式，对原始数据中的异常值进行识别；

所述的异常值检验公式如下：

式中，e_i为修正误差贡献率，θ_i为第i次的体积含水率，为经过第一步至第三步拟合的第i次的预测体积含水率，η_i为第i次的修正系数；

所述的修正系数计算公式为式中，h为压力水头，n为拟合最优的经验参数；

当e_i大于0.1时，表示第i次原始数据存在异常，可以再次试验进行验证。

进一步的，在步骤(1)、(2)及(3)中，所使用的VG模型如下所示：

其中，θ为体积含水率；h为压力水头；θ_r为残余含水率；θ_s为饱和含水率；α和n为经验参数；θ和h为已知参数；θ_r、θ_S、α和n为待定参数。

进一步的，在步骤(1)(2)及(3)中，针对VG模型，使用MATLAB中lsqcurvefit的拟合方法对已知数据进行拟合；具体设置包括初值的设置、待定参数上下限的设置、lsqcurvefit拟合函数的使用，具体的拟合代码如下：

param0＝[0.00001,1.5,0.001,0.3050]；

lb＝[0,1,0,0]；

ub＝[5,5,1,1]；

param＝lsqcurvefit('curvefun',param0,h_data,theta_data,lb,ub,options)

param0为初值设置；中括号内按照前往后的顺序，依次为VG模型中经验参数α、经验参数n、残余含水率θ_r和饱和含水率θ_S；

lb和ub分别为下限设置和上限设置。中括号内按照前往后的顺序，依次为经验参数α、经验参数n、残余含水率θ_r和饱和含水率θ_S；

param为经验参数α、经验参数n、残余含水率θ_r和饱和含水率θ_S的拟合值；

lsqcurvefit为MATLAB中基于最小二乘法的内置拟合方法；

curvefun为待拟合的具体函数，此处函数为VG模型。

进一步的，在步骤(1)(2)及(3)中，所述拟合优度用于度量拟合曲线对于原始数据拟合效果的好坏，其公式如下所示：

其中，R²表示拟合优度，y表示待拟合的数据，表示y的均值，/>表示拟合数据。

实施例：

本发明所述的一种针对VG模型的参数拟合方法，包括以下步骤：

(1)、采用日立CR21N型高速恒温冷冻离心机获取塑性混凝土土样1和土样2在不同转速下的体积含水率，如下表1所示：

表1不同转速下体积含水率

由于塑性混凝土属于混凝土类，选择拟定初值如下表所示；

表2第一步拟合初值

	α	n	残余含水率	饱和含水率
					1	0.00001	1.5	0.001	0.3067
2	0.00001	1.5	0.001	0.3377

选定拟定初值，确定编号1和2土样的饱和含水率初值，根据500rmp、1000rmp、3000rmp、5000rmp、7000rmp、9000rmp和10000rmp下的数据进行拟合，以体积含水率为因变量，水头h为因变量，得到各个参数的最佳拟合值；土样1的饱和含水率最优拟合初值0.3051，拟合优度0.9621；土样2的饱和含水率最优拟合初值0.3350，拟合优度0.9814；

(2)、将上一步至拟合最优的饱和含水率设为初值，将饱和含水率拟合最优时的经验参数n和经验参数α进行调整后设为初值，对已知数据进行拟合，直至经验参数n和经验参数α拟合达到最优；如下表3和表4所示：

表3第二步拟定初值

	饱和含水率	α	n	拟合优度
					1	0.3051	0.000037337	1.6504	0.9621
2	0.3350	0.000070864	1.3095	0.9814

表4第二步拟合最优值

	饱和含水率	α	n	拟合优度
					1	0.3051	0.000033839	1.7971	0.9691
2	0.3376	0.000083848	1.1967	0.9838

(3)、将第二步最终得到的饱和含水率、经验参数n、经验参数α和残余含水率进行调整后设为初值，对已知数据进行拟合，直至饱和含水率、经验参数n、经验参数α和残余含水率拟合至最优，如下表5所示：

表5拟定最优值

经过第一步的初值设置拟合，土样1和土样2的拟合优度分别为0.9621和0.9814；由此可以得到，根据不同土体分类，确定不同的初值的方法效果良好；经过三步的拟合优化，土样1的拟合优度从0.9621增加到0.9691，土样2的拟合优度从0.9814增加到0.9838；从图片2和图片3中可以看到，优化后的曲线更加贴近真实数据值；由此可以得到，经过第二步和第三步的优化，拟合优度在稳步地提高，逐渐接近拟合最优解；

(4)、通过异常值检验公式，计算得到编号1和编号2的不同转速原始数据的误差贡献率，如下表所示：

表6修正误差贡献率表

	500rmp	1000rmp	3000rmp	5000rmp	7000rmp	9000rmp	10000rmp
								土1	0.6312	0.0174	0.0000	0.2267	0.0027	0.0870	0.0350
土2	0.6982	0.0197	0.1331	0.0266	0.1021	0.0202	0.0001

由于VG模型函数本身的非线性，水头高度越大，得到的体积含水率和原始值的误差会越大；单从两者之间的差值去判断原始值的异常，并不合理；故而对误差进行修正校准，提出异常值检验公式。

针对土样1,500rmp和5000rmp的修正误差贡献率为0.6312和0.2267，故而土样1的500rmp和5000rmp存在异常，可以考虑再次试验校准；针对土样2,500rm、3000rmp和7000rmp的修正误差贡献率为0.6982、0.1331和0.1021，故而土样2的500rmp、3000rmp和7000rmp存在异常，可以考虑再次试验校准。

经过四步的优化处理(迭代次数大约6次)，实现了对VG模型的拟合优化，并识别出原始数据中的异常值；而在没有规则的情况下盲猜输入初始值进行优化所需时间较长，且最终拟合值最优的可能性较低。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (9)

1.一种针对VG模型的参数拟合方法，其特征在于，具体实现步骤如下：

(1)、根据不同土的分类，选取VG模型中饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r四个参数的初值，对已知土体的吸力水头h和体积含水率θ的数据进行拟合，从而得到饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r的拟合值，

将初值中的饱和含水率θ_S用得到的饱和含水率θ_S拟合值进行替换，保持其他三个参数的初值不变，后再进行拟合；直至第N+1次的饱和含水率拟合值与第N次的饱和含水率拟合值之差的绝对值小于0.0001，且第N+1次的拟合优度与第N次的拟合优度之差的绝对值小于0.0001；

(2)、将步骤(1)中最终进行拟合得到的饱和含水率θ_S、经验参数n和经验参数α的拟合值替换步骤(1)中选取的三个对应参数的初值，保持残余含水率θ_r的初值不变，对已知数据进行拟合，直至第M+1次的经验参数n的拟合值与第M次的经验参数n的拟合值之差的绝对值小于0.002，直至第M+1次的经验参数α的拟合值与第M次的经验参数α的拟合值之差的绝对值小于0.0000001，且第M+1次的拟合优度与第M次的拟合优度之差的绝对值小于0.00001；

(3)、将步骤(2)最终进行拟合得到的饱和含水率θ_S、经验参数n、经验参数α和残余含水率θ_r的拟合值替换步骤(1)中选取的四个对应参数的初值，对已知数据进行拟合，直至第i+1次的残余含水率θ_r的拟合值与第i次的残余含水率θ_r的拟合值之差的绝对值小于0.0001，且第i+1次的拟合优度与第i次的拟合优度之差的绝对值小于0.000001；

(4)、将步骤(3)最终进行拟合得到的拟合数值和原始值进行对比，通过异常值检验公式，对原始数据中的异常值进行识别。

2.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(1)中，所述不同土的分类共分为四大类，包括普通土体、砂土类、粘土类和混凝土类；

所述普通土体包括粘壤土、壤土、粉砂及粉壤土；

所述砂土类包括壤质砂土和和砂；

所述粘土类包括粉质粘土、粉质粘壤土和砂质粘土；

所述混凝土包括塑性混凝土和混凝土。

3.根据权利要求2所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(1)中，所述土体的饱和含水率θ_S初值根据低转速下的体积含水率确定，取最低转速下的体积含水率为饱和含水率θ_S的初值、且保留四位有效数字；

其中，所述普通土体的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.01、1.4和0.001；

所述砂土类的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.01、2.2和0.001；

所述粘土类的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.01、1.1和0.01；

所述混凝土的残余含水率θ_r、经验参数n和经验参数α的初值分别为0.001、1.5和0.00001。

4.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(2)中，所述的经验参数n保留四位有效数字设置为下一次拟合初值；所述经验参数α保留三位有效数字设置为下一次拟合初值。

5.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(3)中，所述残余含水率θ_r保留一位有效数字设置为下一次拟合初值。

6.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(1)、(2)及(3)中，所使用的VG模型如下所示：

其中，θ为体积含水率；h为压力水头；θ_r为残余含水率；θ_s为饱和含水率；α和n为经验参数；θ和h为已知参数；θ_r、θ_S、α和n为待定参数。

7.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，在步骤(1)(2)及(3)中，针对VG模型，使用MATLAB中lsqcurvefit的拟合方法对已知数据进行拟合。

8.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(1)(2)及(3)中，所述拟合优度用于度量拟合曲线对于原始数据拟合效果的好坏，其公式如下所示：

其中，R²表示拟合优度，y表示待拟合的数据，表示y的均值，/>表示拟合数据。

9.根据权利要求1所述的一种针对VG模型的参数拟合方法，其特征在于，

在步骤(4)中，所述的异常值检验公式如下：

式中，e_i为修正误差贡献率，θ_i为第i次的体积含水率，为经过第一步至第三步拟合的第i次的预测体积含水率，η_i为第i次的修正系数；

所述的修正系数计算公式为式中，h为压力水头，n为拟合最优的经验参数；

当e_i大于0.1时，表示第i次原始数据存在异常，可再次试验进行验证。