土工膜筒状双向拉伸试验非接触式应力应变量测方法

杨 帆 束一鸣 满晓磊 顾 克 郝雪航 徐红霞

（1.河海大学 水利水电学院，南京 210098；2.合肥工业大学 土木与水利工程学院，合肥 230009）

土工膜相比混凝土、粘土等防渗材料，具有防渗效果好，渗透系数仅为10－13m／s，几乎不透水；结构简单、便于施工铺设 ；材料为高分子聚合物柔软材料，单向拉伸中延伸率可达700%甚至更大［1－2］，相比混凝土坝更能适应地基沉降带来的变形.我国的水能资源在西南地区极为丰富，正处于快速发展时期.但是我国西南地区地基多为深覆盖层地基，面板堆石坝在深覆盖层上时混凝土面板容易开裂，防裂工艺不是十分成熟，故而土石坝成为唯一适合在深覆盖层地区筑坝的大坝.土石坝防渗传统多采用粘土防渗，但是粘土开挖对周围生态环境影响较大.土工膜渗透系数非常低、延伸率高的特点使其成为理想的深覆盖地区高坝防渗材料.然而土工膜在大坝运行过程中，承受水压力荷载、地基沉降的影响、温度变化的应力都使土工膜处于双向拉伸应力状态［3］.而关于土工膜的研究成果基本上是基于单轴拉伸，但是在双向拉伸状态下的极限延伸率、屈服强度等力学特性与单向拉伸相比有着显著的降低，以聚乙烯土工膜为例，其纵向断裂强度下降21.53%，横向断裂强度下降60.32%［4］，极限延伸率的变化与加载速率和温度有关，但是在双向拉伸的状态下，极限延伸率相比同条件的单轴拉伸也有显著降低［5］，这种降低显然对于实际工程的设计安全系数影响重大.故而，土工膜双向拉伸试验对于土工膜应用于大坝防渗有着重要的意义.土工膜筒状拉伸试验简单易行，应力分布较平面十字形双向拉伸更为均匀［6］.但是由于圆筒状膜应变较大，应变片是利用导体半导体拉伸过程中的电阻变化来测量金属类小变形的应力测量方式，本次试验需要测得土工膜极限延伸率，变形高达30%以上，通过贴应变片等方式进行测量可靠性不高；且试验材料为立体的圆筒形，难以像平面十字形双向拉伸试验那样进行平面的直接测量.鉴于以上问题，采用数字图像处理技术［7－9］测量各向应变进而计算应力［10－13］，准确测量土工膜在筒状拉伸试验中的变形进而得到其应力应变.

1 筒形双向拉伸试验的设备及原理

根据试验采用的夹具大小裁剪合适的土工膜，将所裁土工膜用进口焊机焊接成筒形，将焊接成的土工膜薄壁圆筒用四条钢箍固定在与单向拉伸万能试验机连接的两个圆盘式的上下夹具上.其中，上夹具上的进气孔通过调压阀与空气压缩机相连，下夹具上的出气孔与精密气压表相连.为实现土工膜的环向拉伸，试验过程中用空气压缩机向圆筒内充气，并通过读取气压表的读数利用调压阀进行反馈调节以保证筒内气压稳定.为了保证圆筒内的密封性，在圆筒与夹具之间的止气槽内填充了橡胶皮条及止气带进行止气.同时由电脑控制万能试验机按一定速率对土工膜进行轴向拉伸，由此实现了土工膜的双向拉伸.具体设备见图1.

图1 筒形双向拉伸试验设备

2 图像采集与处理

数字图像处理技术是指利用计算机、高清相机与各种软件，对图像信息对应的电信号进行相应的运算，例如图像滤波、分割图像、提取特征、图像增强、图像复原等，以便获得需要的特征信息，方便后期运算.为了避免专业基础程序的编写，试验使用Matlab来进行大部分的图像处理.

2.1 图像采集设备

试验图像采集的主要设备为单反相机，如图2所示.为了方便后期进行有效地图像识别，在获取图像时应增强对象与背景的对比度，因此，在薄壁圆筒的后面设置了一片2m×4m的黑幕.为了保证图像位于同一位置拍摄，相机用三角架固定于薄壁圆筒正前方.此外，拍摄过程均使用快门控，避免了直接按动快门键引起相机的微小抖动.拉伸开始前拍摄第一张图像作为初始图像，如图3所示.

图2 图像采集主要设备

图3 初始图像

2.2 图像采集设计

经过多次预实验，可以发现土工膜都是在圆筒中部试验横截面最宽处首先达到屈服，如图4所示.

图4 土工膜屈服状态

经本课题组姜晓桢博士分析，土工膜环向应力的大小沿膜膨胀后的轮廓各处并不相等，而在其最大横截面处最大；轴向应力大小沿土工膜膨胀后的轮廓方向不发生变化.所以，当使用筒状试样作为双向应力条件下的拉伸试样时，最大截面处将首先发生屈服，观测的着眼点也应在最大横截面处，而最大横截面处即为土工膜膨胀后轮廓的最宽处.综上，土工膜筒形双向拉伸试验中应力应变的量测与试样膨胀后的轮廓紧密相关，所以图像采集及处理的重点就是测量试验中试样变形轮廓.

2.3 图像校正

拍摄的图像可能会因设备和人为等复杂原因产生失真.在图像的采集过程中，一般来说图像失真主要来自两个方面：桶形失真和梯形失真.桶形失真是指由成像画面呈桶形膨胀状的现象.在变焦镜头和广角镜头的最广角端，桶形失真比较其明显，这种失真是由于相机自身硬件缺陷引起的.梯形失真则主要是拍摄时镜头与筒状土工膜没有保持平行引起的，是人为因素影响的.为了尽可能消除图像失真对后续数据计算的消极影响，有必要对图像进行复原处理以还原它的真实信息，即图像校正［14－15］.试验采用了 Photoshop Cs6中自带的对EF－S 18－135／3.5－5.6IS镜头的专用镜头校正滤镜（Lens Correction filter），对采集到的图像进行校正，复原效果见图5，其中图（a）、（b）分别为校正前后的图像，由图可以看出，复原效果比较明显.

图5 镜头校正复原效果

2.4 图像增强

由于系统或者环境噪声、曝光不足或者过量等原因，数字图像在其获取和传输过程中会受到干扰，使图像的质量降低，在机器视觉处理中增加了错误率，对于图像的处理极为不利.这时候就需要对图像进行增强处理，以满足人对于图像的视觉需求，增强机器对于图像数据的感知能力.

2.4.1 图像灰度化

先对图像进行灰度化处理，即图像由相机拍照的RGB彩色图像对应的三维矩阵，变为灰度图像对应的二维矩阵，图像失去色彩（即饱和度）会便于数据运算而丝毫不影响对于变形的测量，灰度即图像亮度，灰度值0即为黑色，灰度值为255即为白色，灰度值为0到255之间的为亮度不同的灰色.图像的灰度化选用了Matlab工具箱中的rgb2gray函数来进行.

2.4.2 图像灰度变换增强

为了便于计算机处理后续数据，需要对经过灰度化的图像进行灰度变换加强处理，这样可以提高计算机和人的视觉效果.灰度变换方法很多，最简单有效的就是直方图变换.图像直方图就是图像灰度分布密度函数图，直方图是一个离散的函数图像，横坐标为灰度级，纵坐标为该灰度出现的次数.在灰度直方图中找到分割较暗和较亮的阈值，由图6可知，阈值为50.使用imadjust函数调整对比度，灰度范围［0，255］对应imadjust函数参数为［0，1］，50即为0.2，将函数参数中的low＿in设置为0.2，high＿in为1，low＿out为0，high＿out为1.低于灰度值50的将被设置为灰度值0即为黑色，low＿in到high＿in与增强后的low＿out到high＿out相匹配的关系如图7所示.

图6 灰度直方图

图7 灰度变换前后关系图

使灰度变换前后的图像如图8所示，图8（a）和图8（b）分别为变换前后的图像.

图8 灰度变换增强效果

2.5 提取特征边缘

图像的边缘是指周围像素灰度变化较大的像素点的集合，在对图像进行提取边缘的过程中，经常会受到噪声的不利影响，如果直接对图像提取边缘，常会把噪声误认为是边缘提取出来.因此要求采用的算法能够较好地消除噪声的不利影响，去掉图像的伪边缘，准确识别图像的边缘信息.试验分别采用了Sobel、Canny、Roberta、Laplacian of Gaussian 4种算法进行了边缘提取作为对比，最后确定采用Canny边界检测算子并综合数学形态学中的腐蚀及膨胀等运算进行边缘提取，达到了预期效果，如图9所示.

图9 边缘

3 应力应变计算

筒状双向拉伸的轴向应力主要是来源于轴向外荷载作用；环向应力来源于内部气压的扩张作用.由材料力学可知最大横截面处的环向和轴向应力的计算公式分别为

式中，σa为轴向应力，σr为环向应力，P为内部气压，R为最大横截面半径，δ为土工膜的厚度，F为万能拉伸机测得的轴向外荷载.而最大横截面出的轴向应变和环向应变的计算公式分别为

式中，εa为轴向应变，εr为环向应变，l为拉伸前上下夹具间土工膜的长度，L为拉伸过程中上下夹具间变形后的土工膜长度，r为拉伸前圆筒的最大横截面半径，R为拉伸过程中圆筒的最大横截面半径.受力图如图10所示.

3.1 初始准备工作和程序简介

图10 筒状土工膜双向拉伸后形态及受力图

由于夹具是固定的，r自然是不变的已知数据，l可以通过直尺直接测量.需要注意的是，Matlab以矩阵对应图像的像素，但是像素和实际尺寸不是直接对应的.在得到初始值r和l后，需要在第一张图片中，读取图像中上下夹具之间行坐标之差l′和圆筒的直径2r′，得到l，和l、2r′和2r之间的比例关系，本组试验比例系数为0.3.为方便起见，下文中的以l和r均以图像像素值计算，换算成实际长度只需在程序结尾乘以本次实验的比例系数即可.由于每组试验位置不变，此比例系数也是不变的，适用于整组试验.

本次测量轴向和环向的应力应变是基于Matlab编写脚本文件实现的，现简略介绍此脚本文件中编写的代码结构及其功能如图11所示.

3.2 环向应变计算

环向应变的计算关键是提取边缘图像图9中的最膨胀处的列坐标.由于土工膜是可以近似各向同性的材料，同时由实验拍摄的照片可以印证，图像近似于对称图像，为了计算机运算方便，提高运算效率，可以统一以边缘图像左半侧部分为研究对象.将边缘图像中灰度值为1的点通过Matlab编写脚本程序赋值到一个新的矩阵A中，［RC］＝A，即行坐标赋值给R，列坐标赋值给C，如本图中灰度值为1的点有2 888个，R和C即为2 888乘1大小的矩阵，A即为2 888乘2大小的矩阵，A存储了所有边缘点的坐标值.在Matlab中一幅图像对应一个矩阵，矩阵的行坐标即为图像的纵坐标，矩阵的列坐标即为图像的横坐标，故而Matlab以图像的左上角为像素坐标的起始值（0，0）.此时搜索最膨胀处点的纵坐标即为搜索矩阵A中列最小值Cmin.第一步是通过ginput读取第一张图中上夹具底端与土工膜夹持边界的两个端点处的列坐标C1和C2，如图12所示，图（a）中所读取的点为C1，图（b）中读取的店为C2.因为相机在同一组实验中位置不变且万能试验机的上夹具不同，移动下夹具以拉伸材料故而可保证同一组试验中拍摄的图片上夹具的位置（即图像中的坐标不变），C1和C2的中点也保持不动，故而可以通过读取一张图片的C1，C2写入脚本中后就可以实现对于一组实验图像的批处理.第二步就是搜索每张图片中C矩阵中的最小值Cmin，通过min函数实现.每组实验中的R值即为

代入公式（4）即可得到环向应变.

图12 读取C1和C2的坐标图

3.3 轴向应变计算

3.3.1 有效区域分割

轴向应变计算的关键是得到实验过程中拍摄的每张图像的L值.同样地，基于对称性和计算方便考虑，在去掉固定不动的上夹具部分后，通过搜索Cmin得到对应的横纵坐标后，以上夹具底端对应的行坐标即C1和C2对应的横坐标、Cmin对应横坐标截取横向，以纵坐标0和C1、C2中点处纵坐标值截取纵坐标.截取的图像即为有效区域的1／4为计算对象.通过分割边缘图像图9得到图13.

图13 1／4图

3.3.2 变形后曲线长度的计算

由于相机和Matlab以矩阵各点存储图像各点的特征信息，本身即是以矩阵的离散点描述实际图像，而在提取边缘的过程中又依据只提取灰度突变的点作为边缘的原理，故而即使使用高清相机拍摄图像，在相机本身缺陷和Matlab提取边缘的算法原理下，势必出现边界曲线在放大的情况下不够光滑的现象.曲线实际在细部抖动而非平滑，故而应当采取措施平滑曲线，减小应力应变计算的误差.先在Matlab中编写脚本程序M文件，提取图像13中所有灰度值为1的点，行坐标存为R，列坐标存为C，通过差分函数diff依次对R和C求其差分，并平方后求和再开方运算，再求和可得到图13的白色线段之长.现在为求平滑曲线的L／2，每隔5个点取为新的差分形式，即编写循环语句对R和C每隔5个点抽出进入新的矩阵Rxin和Cxin，对Rxin和Cxin进行上述求曲线长度的运算，求取新的L／2长，新L／2会以每6个相邻点的直线长度求和代替每两个相邻点直线长度求和来得到总的曲线长度.使用plot函数分别描绘［RC］和［RxinCxin］即平滑前后曲线，如图14所示，可见曲线变得平滑.以本图为例，图14的曲线（a）长度为212.87即图13中曲线长度，图14（b）的的曲线长度为205.89，使用软绳实测为204，可见平滑方法更接近真实值，且误差不到1%.

图14 曲线平滑前后效果

同样，在一组试验中，由于相机位置不变，截取1／4图的4条直线中的上横线不变，下横线可以通过编写程序搜索列最小值点的行坐标得到，而左边纵直线为定值0，右边纵直线为C1、C2的中点的列坐标，也是定值.故而可以批量处理一组试验中所有图像的L值，然后代入公式（3）中即可求得轴向应变.

3.4 轴向和环向应力的计算

与万能试验机相连的电脑通过拉力传感器可以实时测量每个时刻的拉力F，而环向气压可以通过连接存储型数字气压计记录每个时刻的内部气压，通过USB接口可以很方便地导入电脑中，以表格形式存储.Matlab可以读入EXCEL中的气压值和拉力值，根据之前通过图像处理分析得到的轴向L和环向变形R，而厚度δ在各个时间点的数值由数显测厚规测量，实时计入EXCEL表格中.由公式（1）和公式（2）可以十分便捷地得到环向应力和轴向应力.由此，可以实时测量筒状拉伸试验中的轴向环向应力应变.在40kPa和50kPa两个不同环向压力的情况下，轴向应力应变关系如图15所示.

图15 轴向应力－轴向应变曲线

由于轴向拉伸速率一定，轴向应力随着轴向应变的增大而增大，并且增大的速率逐渐减小，曲线慢慢趋于水平，最终轴向应力不再发生较大变化.比较两张图可以看出，在不同内压下轴向应力所能达到的最终值并不相同，并随着气压的增大而减小.这是因为在双向受拉的情况下，土工膜在相对较低的轴向应力下即达到屈服.

4 结 语

本文进行了不同内部压力下土工膜内压薄壁圆筒试验，运用数字图像处理技术，应用Matlab软件强大的图像处理功能，对土工膜筒形双向拉伸试验过程中尺寸较大且为非平面的变形进行了精确测量，进而得出了不同环向应力条件下的土工膜筒状试样的轴向应力与轴向应变关系曲线及环向应变与时间的关系曲线，通过对试验数据的分析，得到了土工膜在双向应力条件下的变形特点.而所测算得到的各向应力及应变也为进一步探究土工膜本构关系提供了准确的数据.不足之处在于，由于相机像素太高，Matlab运算会变慢，而像素偏低又会出现曲线不光滑的现象，虽然本文采用分段直线差分的方法来平滑曲线，也能达到实验数据要求的精度，但是此处仍然有着改进的空间.

［1］ 张思云，鲁伟涛，靳向煜.PE土工膜双向拉伸试验研究［J］.隧道建设，2012，32（5）：654－657.

［2］ 林继镛，王光纶.水工建筑物［M］.北京：中国水利水电出版社，2009.

［3］ 吴海民，束一鸣，曹明杰，等.土工合成材料双向拉伸多功能试验机的研制及初步应用［J］.岩土工程学报，2014，36（1）：170－175.

［4］ 张思云，张 艳，靳向煜.土工膜和非织造土工布单向与双向拉伸机理对比试验研究［J］.东华大学学报：自然科学版，2014（2）：016.

［5］ Bray J D，Merry S M.A Comparison of the Response of Geosynthetics in the Multi－axial and Uniaxial Test Devices［J］.Geosynthetics international，1999，6（1）：19－40.

［6］ 任泽栋，满晓磊.堆石坝防渗土工膜薄壁圆筒双向拉伸性能测试方法［J］.三峡大学学报：自然科学版，2013，35（3）：21－25.

［7］ 王耀南，李树涛.计算机图像与识别技术［M］.北京：高等教育出版社，2001.

［8］ 余成波.数字图像处理及 MATLAB实现［M］.重庆：重庆大学出版社，2003.

［9］ 张博夫，梁凯琦.MATLAB环境下的数字图像处理实验入门［J］.实验科学与技术，2008，6（1）：52－53.

［10］Fernlund J M R.3－D Image Analysis Size and Shape Method Applied to the Evaluation of the Los Angeles Test［J］.Engineering geology，2005，77（1－2）：57－67.

［11］Chetana B.Rao.Development of Three－dimensional Image Analysis Techniques to Determine Shape and Size Properties of Coarse Aggregate：［dissertation］.Urbaba－Chanpaign：Univ.of Illinois，2001.

［12］张婧娜.基于数字图像处理技术的沥青混合料微观结构分析方法研究［D］.上海：同济大学，2000.

［13］Ramitha Wettimuny，Dayakar Penumadu.Application of Fourier Analysis to Digital Imaging for Particle Shape Analysis［J］.Journal of Computing in Civil Engineering，2004，18（1）：2－9.

［14］姜大志.标准图形法求解相机镜头非线性畸变的研究［J］.东南大学学报：自然科学版，2001，31（14）：111－116.

［15］李 嘉，林 辉.基于数字图像处理的粗集料棱角性量化研究［J］.公路交通科技，2008，25（7）：27－31.