杜瑞锋, 裴向军, 张晓超, 常志璐, 崔雪婷

(1. 成都理工大学 地质灾害防治与地质环境保护国家重点实验室， 成都 610059；2. 内蒙古建筑职业技术学院 建筑工程学院， 呼和浩特 010070)

黄土击实试验曲线峰值点的数值方法

杜瑞锋1,2, 裴向军1, 张晓超1, 常志璐1, 崔雪婷1

(1. 成都理工大学 地质灾害防治与地质环境保护国家重点实验室， 成都 610059；2. 内蒙古建筑职业技术学院 建筑工程学院， 呼和浩特 010070)

击实试验是衡量实际施工中土体被压实或夯实到最优程度的室内指导性试验，对土木工程的抗压、抗剪、抗渗等性能有重要意义。在探讨击实试验中的密度、含水量、干密度在数学及土力学意义的空间图形关系基础上，应用数值分析原理并采用分段三次多项式插值方法，自编Matlab程序，有效、快捷地确定击实曲线峰值点以及绘制成果曲线，克服了应用“图解法”的不准确性和不方便性。该方法可应用于不同地区的黄土击实试验中，效果良好。

击实试验； 图解法； 峰值点； 数值方法

0 引 言

击实试验是土木工程中的基本土工试验之一，试验目的在于确定土体在某一击实功下的最优含水率和最大干密度，用于评价建筑地基土、路基土、挡土墙回填土等的压实或夯实质量，对土体的抗压、抗剪、抗渗等性能有重要意义[1-2]。因此，击实试验结果准确性对工程建设方、施工方、监理方等单位都是重要的、敏感性的指标。

多数研究者从击实试验结果的影响因素、如何提高试验结果的精确度以及不同类型土的击实特性等方面进行探讨[3-4]；极少数研究者对击实试验数据处理问题进行过研究，但缺乏系统的、科学的、实用的研究方法和工具[5]。

黄英等[6]于2002年提出通过数值解法获取击实试验基本参数，认为最大干密度和最优含水率之间符合2次函数关系，对试验数据进行拟合处理，可用于数据点不足5点的情况下的分析，由于采用了数值分析中的拟合处理，误差是不言而喻的，实际中应用起来不方便，且没有给出实用的工具。

王永和等[7]于2006年利用Matlab软件对击实数据进行了拟合处理，分别采用了2次、3次多项式进行拟合，认为3次多项式曲线能较好地拟合最大干密度和最优含水率之间关系，同样采用了拟合方法，但缺乏有力的分析和实用的程序工具。

本文针对黄土的压实特性研究中曾做过多组轻型击实试验。在整理试验数据时，常遇到曲线不理想而需要补点或重做现象；而更关键的问题是如何科学、统一、方便地确定击实曲线的峰值点。因此，如何使击实试验数据处理过程既简单又不失科学性，是本文探讨的主题，即寻求一种既省时、省力，又能快速、准确地确定土的最优含水率和最大干密度的数据处理方法。

1 击实试验

击实试验目的在于确定某种土的最优含水率和最大干密度。当击实功和击实方法一定时，土的干密度随含水率的增加而增大；当干密度达到某一最大值后，含水率的继续增加反而会使干密度减小。因此，干密度与含水率之间存在着一定的函数关系，该函数的峰值点即被确定为土的最优含水率和最大干密度。

击实曲线之所以呈现曲线变化的特点，是因为土中的水分所处的状态。当土处于低含水率时，土中水分主要是结构水或强结合水。击实时，由于颗粒表面基本无薄膜水，随着干密度增大，摩阻力增大，击实功较多消耗在克服土体移动的摩阻力上，有效击实功降低，因而土体不易被压实。

当含水率逐渐增大时，颗粒表面水膜逐渐变厚，会出现大量的薄膜水，其水膜的润滑作用增大。击实时，颗粒表面摩阻力相应减小，有效击实功增加，压实效果增强。但当含水率超过某一值后，土颗粒中孔隙内的自由水逐渐增多，大量自由水被封闭不易排出，击实功将大量消耗在自由水中，击实效果又将降低，同时大量孔隙被自由水占据，土的干密度就会减小[2，8]。

2 击实试验数据

2.1 击实试验数据处理方法

我国不同行业的土工试验规范或规程中对击实试验结果的处理原则是以干密度为纵坐标、含水量为横坐标，曲线上峰值点的纵、横坐标分别为最大干密度和最优含水量[9-11]。这里首先需要说明的是，曲线上的峰值点不能简单地理解为是击实数据中的最高点，注意强调的是曲线上的峰值点。

在过去没有计算机普遍应用的年代，击实试验数据处理一般是在坐标纸上手工描绘出至少5个点的数据，然后用光滑曲线连接起来，然后找出曲线上的峰值点，定为最大干密度和最优含水量，有文献称之为图解法。显而易见，局限于手工绘制光滑曲线的缘故，最终确定的曲线上峰值点常常就是坐标纸上描出的最高点。这就使得击实试验数据处理过程中结果随意性大，存在人为的误差等等，失去了科学的数据处理原则，同时忽略了数据变量间的内在的数学函数关系。

如图1(a)所示，表面上看，6个点的试验数据中第5个点的位置最高(w=15.8%，ρd=1.83 g/cm3)；如果采用手工处理方法找峰值点，第5点会被大多数土工试验员选中；而图1(a)中给出的峰值点坐标值为：wop=15.0%，ρd,max=1.83 g/cm3，显然这两组数值是不相等的。此处，姑且不论文献[9]中采用了何种方法确定出了曲线上的峰值点，但绝对不是手工“图解法”所能完成的。更进一步推断，图1(a)中击实试验数据处理中使用了数值分析中的插值运算，即对击实数据进行了插值处理。当然，这是基于目前计算机普遍应用的前提下能完成的工作，对土工试验员或工程师而言，当然也需要一定的数值分析基础知识。顺便指出，击实试验成果曲线中，应绘出饱和曲线(图1(a)或(b)中右侧的倾斜直线)，并标以Va=0，即表示土体中空气体积为零。

(a)文献[9]中曲线(b)本文方法绘制的曲线(o为峰值点)

图1 击实试验曲线

很显然，上述“图解法”不能满足当今科研和生产实践中的要求。为了适应计算机处理数据的需要，对击实数据整理可在数值分析原理的基础上编写适用的计算机程序，以减轻“图解法”带来的繁琐工作，提高数据处理的科学性和准确性。

2.2 击实试验数据的空间关系分析

在土力学中，土的密度ρ、含水量w、干密度ρd之间有下式关系：

如果单纯地在数学空间内表达三者的关系，即如图2(a)所示的曲面，显然函数ρd是不存在极限值的；而在击实试验中，ρ、w是2个独立测量的物理变量，ρd是由ρ、w在土力学理论框架下导出的物理变量，赋予了土力学实际概念且基于文献[9]中数据绘出的3个变量关系散点图如图2(b)所示，如果把ρ、w看成2个连续的变量，那么组成的曲线是有极值的，击实试验就是确定ρ、w的极值问题的试验。将图2(a)、(b)表示的图形绘制在同一坐标下，即为图2(c)所示的图形，可见散点落在曲面上。

(a) 数学空间分布

(b) 土力学意义分布

(c) 散点图与曲面叠加图

图2 土的密度ρ、含水量w、干密度ρd之间空间关系

3 应用Matlab处理击实试验数据

3.1 数值分析中的插值与拟合

科学试验及生产实践中，面对一组数据(x1,y1),(x2,y2),…,(xn,yn)，而函数表达式y=f(x)未知，要从某函数类(如多项式函数、指数函数等)中求得一个函数y′=f(x′)作为y=f(x)的近似函数，这类计算问题统称为数据建模或数值分析[12]。

处理试验数据时采用的数据分析方法有两类：① 插值方法，要求所求的函数严格通过数据点；② 拟合方法，函数可以不全通过数据点，即允许存在误差，但要求达到某种误差指标的最小化。通常，插值方法比较适合数据准确或数据量较小的情况；而拟合方法适用于数据有误差且数据量大的情况[13]。

击实试验中，由于数据少，一般仅有5、6个数据点。因此，对击实试验数据进行处理，合适的方法是插值方法，而不是拟合方法，即达到使数据点间的曲线符合光滑、平顺的目的。在处理本研究项目中不同地点黄土击实试验时，本文采用分段的3次多项式函数插值方法，可保证每个数据点处(亦称“断点”)具有连续的2阶导数，这样可使得插值函数曲线具有比较光滑、平顺的特性。

3.2 处理程序在Matlab软件平台上实现

基于数值分析原理，在Matlab软件平台上，将不同的击实试验数据写在扩展名为.m的脚本文件中，通过执行自编的M函数程序SKLGP_CompactionTest.m，能完成数据的插值运算、峰值点的确定以及曲线的绘制工作。为了便于交流、讨论及应用，列出了所编写的脚本文件和M函数程序的全部代码。限于篇幅，M程序中没有列出分段三次多项式的方程式，有兴趣的读者可自行处理；同时说明的是，数值分析计算原理在其他文献中有充分的证明过程，因此本文不作相关的讨论和证明。

脚本文件Gonglulizi.m的内容如下：

clc; clear all; format bank

%以下击实试验数据来自文献[9]

rop=[1.71 1.75 1.80 1.83 1.76];

w=[10.2 11.8 13.0 15.8 19.0];

%Name为工程项目名称

Name='文献[9]中击实试验曲线';

SKLGP_CompactionTest(w,rop,Name)

总体上，M程序中插值方法采用分段三次多项式插值，其中spline为三次样条曲线插值算法，pchip为三次Hermite插值算法，cubic为双三次插值算法[13]。3种插值计算方法有所差异，插值效果自然不同，但均具有较好的光滑性。使用时，在interp1命令选项中设定即可。

SKLGP_CompactionTest.m的内容如下：

function SKLGP_CompactionTest(w,rop,Name)

W=linspace(min(w),max(w),100);

Rop=interp1(w,rop,W,′spline′);

Rop_max=max(Rop);

[y,J]=max(Rop);

plot(w,rop,′o′,W,Rop,′-r′,′linewidth′,6)

grid on

hold on

title(Name,′fontsize′,30,′fontweight′,′bold′)

xlabel(′w(%)′,′fontsize′,25)

ylabel(′ ho_d(g/cm^3)′,′fontsize′,25)

set(gca,′FontSize′,20)

wmax=(1./rop-1/2.7)*100;

plot(wmax,rop,′-m′,′linewidth′,4)

set(gca,′Xtick′,[round(min(w))-1:1:round(max(wmax))+1])

set(gca,′Ytick′,[round(min(Rop))-1:0.05:round(max(Rop))+1])

gtext([′最大干密度′,′=′,num2str(Rop_max),′(g/cm^3)′],′fontsize′,20)

gtext([′最优含水量′,′=′,num2str(W(J)),′(%)′],′fontsize′,20)

plot(W(J),Rop_max,′o′,′Markersize′,20)

gtext([′V_a′,′=′,′0′],′fontsize′,20)

hold off

3.3 不同地区黄土的轻型击实试验数据处理

来自4个地区的黄土轻型击实试验中，采用干法制土，每层击数为27。击实曲线见图3，其中图3(a)、(b)中采用spline插值，图3(c)、(d)中采用cubic插值。其中，图3(a)、(c)、(d)中经本文方法插值处理后峰值点位置就是试验散点中的最高点，能与采用“图解法”确定峰值点保持一致性；而图3(b)经插值处理后出现的峰值点位置显然不是试验散点中的最高点。这是因为，峰值点选取了分段三次多项式插值函数曲线的最大值点，当然从数学意义上舍弃了试验点中的最高点。整个过程符合数值计算原理，曲线光滑、平缓，完全符合文献[9-10]中规定的确定峰值点的基本原则。同时需要说明的是，如果利用“图解法”确定峰值点则非常困难，且不易做到数据统一、操作简单。因此，摒弃“图解法”势在必行。

(a) 河南灵宝

(b) 陕西富县

(c) 甘肃兰州

(d) 宁夏石碑塬

图3 不同地区黄土的击实试验成果曲线(o为曲线峰值点)

最后，作为一个有力的例证，应用本文方法处理文献[9]中的击实试验数据。击实曲线如图1(b)所示，与图1(a)相比，尽管两图的纵横坐标比例有差异，但曲线的形态、峰值点位置及数值能保持一致性，表明自编程序的科学性、有效性和可操作性。

4 结 论

(1) 本文指出“图解法”在处理击实试验数据时的不足之处。在目前计算机普遍应用的情况下，科研和生产实践中应摈弃“图解法”；建议在相关的击实试验规范或规程中应明确指出，并建议采用相应的插值方法等；从相关的外文文献中，也可见数值方法确定峰值点的应用实践[15]。

(2) 为满足黄土击实试验的需要，基于数值分析原理和Matlab软件平台，编写了能进行分段三次多项式曲线插值的程序，能科学地、方便地确定出峰值点以及绘制击实成果曲线。为了便于土木工程师、高校师生及土工试验员进一步交流、讨论及应用本文中的自编程序，列出完整的程序代码，以期得到更好的完善。

(3) 黄土是一种大孔隙、结构性强的粉土，有时成分接近于粉砂，但性质又不同于一般的粉土[16]。本研究中仅进行了黄土的击实试验数据处理工作，对于其他类型的土(如黏土以及粗粒土等)是否适用本文提出的插值方法及程序，还需要进一步的试验和探索。

[1] 陈 辉，雷胜友. 土的击实试验影响因素[J]. 建筑科学与工程学报，2005，22(2)：69-71.

[2] 李广信. 高等土力学[M]. 北京：清华大学出版社，2005:197-200.

[3] 杨小红. 击实试验过程中常见问题分析[J]. 建筑设计，2016，45(9)：31-34.

[4] 王 旭,赵秀绍,艾成刚,等. 提高室内击实试验精确度的方法研究[J]. 华东交通大学学报,2013, 30(5):18-23.

[5] Meehan C L, Cacciola D V. Assessing soil compaction using continuous compaction control and location-specific in situ tests[J]. Automation in Construction, 2017，73：31-44.

[6] 黄 英，符必昌. 确定土的最大干密度和最优含水率的数解法[J]. 岩土工程学报，2002,24(4)：538-540.

[7] 王永和，胡 萍，卿启湘. 影响土样击实试验结果的试验分析[J]. 铁道科学与工程学报，2006，3(6)：31-34.

[8] 朱崇辉，王增红，希罗科夫 B H. 单点击实法击实试验研究[J]. 岩土力学，2012，33(1)：60-64.

[9] JTG E40-2007. 公路土工试验规程[S].

[10] GB/T 50123-1999，土工试验方法标准[S].

[11] 林宗元. 岩土工程试验监测手册[M]. 北京：中国建筑工业出版社，2005:80-84.

[12] 李庆扬，王能超，易大义. 数值分析[M]. 5版. 北京：清华大学出版社，2013：35-46.

[13] 姜健飞，胡良剑，唐 俭. 数值分析及其Matlab实验[M]. 上海：科学出版社，2004：73-104.

[14] 刘国良，杨成慧. Matlab程序设计基础教程[M]. 西安：西安电子科技大学出版社，2012：218-246.

[15] Wong L S, Mousavi S, Sobhani S,etal. Comparative measurement of compaction impact of clay stabilized with cement, peat ash and silica sand[J]. Measurement, 2016, 94：498-504.

[16] 邢国起，肖洪天，李大勇. 潍坊中南部平原地区粉土强度特性试验研究[J]. 南水北调与水利科技，2014，12(2)：62-66.

Numerical Method on Peak Point in Loess Compaction Test

DU Ruifeng1,2, PEI Xiangjun1, ZHANG Xiaochao1, CHANG Zhilu1, CUI Xueting1

(1. State Key Laboratory of Geo-hazards Prevention and Geo-environment Protection，Chengdu University of Technology，Chengdu 610059, China； 2. Department of Civil Engineering, Inner Mongolia Technical College of Construction, Hohhot 010070，China)

In actual construction, compaction test is used to evaluate the tamping quality of soil and reach the optimal state, and is a laboratory experiment with significance to soil performance of compression, shear and penetration. The relations of spatial drawing and bulk density, water content and dry density in compaction test of civil engineering are researched respectively with viewpoints of mathematics and soil physics. A self-made program used Matlab is installed successfully on the platform, and it can fast deal with peak points and plot curve. The method gets better result than the diagrammatical method which is still used now, and is inaccuracy and inconvenience. The program is suitable to other soils (such as clay and coarse-grained soil), and is needed to develop more related research.

compaction test; diagrammatical method; peak point; numerical method

2016-12-20

国家重点基础研究发展计划项目(2014CB744703)

杜瑞锋(1978-)，男，山西五寨人，博士生，主要研究方向为地质工程及土木工程。

Tel.：18686089280; E-mail：dddrrrfeng@126.com

裴向军(1970-)，男，内蒙古林东人，教授，博士生导师。

Tel.：028-84079076； E-mail：peixj0119@tom.com

TU 411.99

A

1006-7167(2017)08-0044-04