光弹性颗粒的线扫描测力方法

陈福星,杨朝宇,邵荣生,黄旋格,张兴刚

(贵州大学 物理学院，贵州 贵阳 550025)

光弹性效应[1]在研究颗粒物质也具有很大的意义. 采用光弹性研究二维颗粒，可通过拍摄光弹性颗粒上的全场光强分布，从而获取颗粒在偏振光场中的干涉条纹分布图，即应力光图. 应力光图反映了二维颗粒的应力分布. 通过应力光图反应的应力分布，可能反推出颗粒边界上的接触力大小. 数字图像的处理[2-3]与二维颗粒的光弹性理论相结合，通过应力光图获知加载力大小及颗粒接触点位置等信息. 光弹性理论的难点在于光弹性模型的应力分布的求解[4]. 对于两点对心对称加载颗粒[5-6]的理论已经有较为全面的研究. 而要处理颗粒的多点不对称接触的复杂情况，则需要更加完善的多点接触圆盘理论[6-7].

杜克大学Behringer等，用光弹颗粒实验探究了力网、力的概率分布、颗粒物质的堵塞态和堵塞转变等问题[8-9]，包括使用对颗粒的应力光图的等差条纹图进行非线性最小二乘法拟合获取颗粒的加载力大小的方法，但该方法是对每个颗粒的全场光强进行拟合，在应用中有较大的运算量，且只有当3维圆柱颗粒的接触方式接近线接触时才较为实用，其对实验拍摄的应力光图的精度要求也较高. 基于平均平方灰度梯度获取平均接触力[10-11]的经验方法，在处理大量颗粒平均接触力时较为适用，但它不能得到颗粒上各个加载点的精确的加载力大小，只能获得颗粒接触力的平均值. 通过对颗粒进行阈值切割从而实现几何形态学操作的方法可以获取颗粒圆心位置、大小、接触点[11-12]. 然而这种方法检测接触点时，存在拍摄角度、图像识别精度等问题，可能会把2个靠得很近却不接触的颗粒识别为相互接触. 这样在计算平均接触力时，由于接触点的数目的误判，会带来平均接触力大小的判断错误. 另外对于单个颗粒，或者颗粒堆边缘处的颗粒，要通过该方法获取接触点位置也是不可能的.

为了解决颗粒阈值切割寻找接触点的位置出现误判，并能寻找边缘处颗粒或者单个颗粒的接触点位置，需要一种新的寻找接触点方法. 而为了解决对颗粒的全场光强进行的非线性拟合方法运行慢、依赖线接触方式的问题，考虑对颗粒全场光强的拟合转换为对颗粒上几条连线的光强进行拟合，即对颗粒的应力光图进行线扫描.

1 二维多点加载的颗粒的光弹性原理

光弹性实验光路如图1所示，在正交圆偏光场中，对于半径为R的二维光弹性材料的颗粒圆盘，存在厚度的三维圆柱，在其圆边界处若存在点加载的载荷，由于光弹性材料的暂时双折射效应，会出现明暗分布的等差条纹. 条纹的光强[1]为

(1)

其中，Ia为光源光强，δ为某处光弹性材料的相位，N为条纹级数，对于整数级的条纹，如N=0，1，2，…，光强I为0. 又比如当I=Ia/2时对应的条纹级数为N=1/4，3/4，5/4，…

光弹性颗粒上每个位置的相位都由其应力分布得到，即二维应力光学定理：

δ=2π(σ1-σ2)h/f,

(2)

其中(σ1-σ2)，h和f分别为主应力差、模型厚度和材料条纹值.

图1 光弹性实验光路图

要求得光弹性颗粒光强的分布，还需要求出多点加载时颗粒圆盘各处的主应力差(σ1-σ2). 求出单点载荷在圆盘上的应力分布，通过应力的叠加原理，即可得到多点加载的圆盘的应力分布.

如图2所示，对于过颗粒模型圆心的加载力F，要求它在圆盘上产生的应力分布，可以建立以加载点为原心的坐标系O1，以及在圆心处建立的直角坐标系O.O1坐标系y轴与O坐标系的y轴夹角为α，在圆心坐标O中，圆盘上任意一点坐标可表示为(r,θ)，r为P点到圆心距离，θ为P点到原心连线与y轴的夹角.

图2 与O坐标系的y轴夹角为α的力作用在颗粒模型上的示意图

同时，P点也可在O1坐标系下的坐标可表示为(r1,θ1)，在O1坐标系中可给出力F在y轴的正半平面中的应力函数[7]为

(3)

为了方便应力叠加，将应力函数在O坐标系中表示. 坐标系O的可由坐标系O1先通过平移R距离再逆时针旋转α变换而来，由几何关系

(4)

则坐标系O下的F产生的应力函数写为

(5)

应力张量可以由应力函数[4]表达为

(6)

如图3所示，现在考虑有i个点加载力作用在圆盘上，若考虑在颗粒圆盘边界处的加载均为无摩擦的对心点加载，现考虑这i个点加载力半平面的应力叠加近似为圆盘上的应力分布，其表达式为

(7)

由于主应力差为

(8)

即可由此得到颗粒圆盘的光强分布. 虽然是近似结果，在模拟中发现其完全能够满足实验的要求.

图3 多个力(F1,F2,F3,…,Fi)作用在颗粒上的示意图

2 加载实验装置

为搭建多点加载颗粒的实验装置如图4所示，可先设计所需的具有不同加载点位置的加载装置的模型，并使用3D打印机打印制作，安装螺旋微分头. 为了能够实现在颗粒上进行多点加载，同时也能够检测出加载力的大小，需要在螺旋微分头与颗粒圆盘之间放置力传感器. 要使力传感器的检测头能与螺旋微分头之间不出现相对滑动，还需要使用3D打印制作小的力传感盒来固定测微头与应力检测头的相对位置. 搭建正交圆偏光场，需要打印能搭载偏振片和1/4玻片的装置，该装置应能实现玻片与偏振片的在一定角度范围内的旋转，以方便调节得到正交圆偏光场. 另外由于显示屏上已有偏振膜，所以不用再添加起偏镜，仅需要旋转调节1/4玻片与检偏镜即可得到正交圆偏光场.

图4 实验装置搭建与加载装置示意图

制作光弹颗粒时，可先使用3D打印机打印凹形圆柱硅胶模具的模具，由于打印精度的问题，在模具圆面会有一定程度的粗糙，可使用细砂纸将其打磨平整. 使用硅胶浇筑，得到光弹性材料的硅胶模具. 另外，为了制作能在力传感器的检测范围内有丰富图案变化的光弹性颗粒，可调试需要的比例，将环氧树脂的软胶和硬胶按需求进行混合，可浇筑得到不同材料条纹值的颗粒试件.

3 对颗粒应力光图的数字图像处理

3.1 颗粒圆盘的边界、圆心、半径的提取

获得颗粒实验的应力光图如图5(a)所示，可以发现，由于二维颗粒存在一定厚度，对环境光源有折反射，在颗粒边界处，也有一定的光强. 可通过简单的阈值二值化[图5(b)]、闭操作以及填充操作获得颗粒圆盘的范围[图5(c)].

(a)实验颗粒应力光图

(b)对实验应力光图进行阈值二值化

(c)闭操作及填充操作得到颗粒形态矩阵图5 颗粒圆盘的边界、圆心和半径的提取

对颗粒圆盘范围图5(c)进行遍历，可以确定颗粒圆盘范围内的像素点，通过质心公式得到圆盘圆心. 统计模型范围内像素点数，即可用圆面积公式得到圆盘的半径[12]. 如此将颗粒圆盘的边界抽象为圆.

3.2 边界圆模板寻找加载点位置的方法

当处理多颗粒接触问题时，可直接采用阈值分割的方法对颗粒几何形态学操作获取颗粒的边界，对2个颗粒边界相切的位置即可直接判断为接触点位置[11-12]. 然而这种方法检测接触点时，存在拍摄角度、光照、图像识别精度等问题，会把2个靠得很近却不接触的颗粒识别为相互接触. 另外有一种情况，虽然颗粒接触但却没有接触力存在，这样把它判断为接触也是不合适的.

本文提出直接通过边界圆模板寻找加载点位置的方法. 直观观测得知，接触点的位置即条纹产生的位置，高级条纹在接触点附近密集，于是其位置对应的像素灰度梯度方就较大，相反地，其他区域的灰度梯度方较小. 如图6(a)，图片某一像素点Ii,j的灰度梯度方[10-11]可由其8个邻域像素点及其自身的光强定义为

(9)

(a)计算Ii,j像素点的平方灰度梯度

(b)实验的灰度梯度方分布图

(c)模板内H与θ的关系图像图6 边界圆模板寻找加载点位置的方法

在处理接触点位置时，发现单单使用边界圆模板法得到接触点位置有时会有1°～4°的不确定度，而接触位置的误差可能会影响最后处理得到的加载力大小，所以，若颗粒数较多时，采用边界圆模板法与阈值切割法[11-12]相结合的方法，先采用边界圆模板法获得加载点的初步位置，再对应阈值分割法得到的精确位置，如此可结合2种方法的优势，既不出现对接触点的误判，也能保证接触点位置的精度. 处理单个颗粒或者少量的边缘处的颗粒时，若想获得更精确的加载点位置，则需要让边界圆模板法获得的结果再经过人为检验来确定更精确的加载点位置.

4 线扫描光强拟合确定接触力的方法

通过3.2的方法获取得到加载点位置后, 通过加载圆盘的应力分布的理论可知圆盘内的光强的分布，从而由光强的变化反推出加载力的情况. 然而，对于颗粒的应力光图，若使用该颗粒全场应力光图作拟合，计算量与误差都很大. 接下来只考虑每个加载点到圆心的连线上的光强变化.

4.1 离散坐标系下直线的确定

由于图片是离散的像素点的灰度值矩阵，要得到加载点到圆心2个像素点的连线的灰度变化，先确定哪些像素点属于在离散的坐标系中2个像素点之间的“连线”，“连线”实际上是一系列的像素点. 为得到两像素点之间“连线”，可参考Bresenham算法[13]. 该算法虽然可以得到2个像素点之间的“连线”，但得到的“连线”的“长度”，即找到的像素点的个数，在两像素点所在直线的斜率为非零有限值时并不合适. 因为颗粒为圆形，那么加载点到圆心2个像素点距离的定义应该按照欧氏距离来定义才较为合理，比如，起点像素点坐标为(10,5)，终点像素点坐标为(40,45)时，“连线”像素点个数按照欧氏距离的定义应该为50个或接近50个，然而使用Bresenham算法得到的“连线”的“长度”为30个或40个.

为了解决这个问题，提出了像素点的逐步行进的方法,如图7所示，从起点像素点开始，对当前像素点十字4邻域的像素点进行判断，判断满足以下2个标准的像素点则为下一步需要走的像素点：1)该像素点与起终点之间的理论直线的点线距离最短;2)该像素点与终点的直线距离最短. 注意，如在像素点的十字4邻域中，某一点为已经走过的点，则需要在判断时将其略去，这样逐步行进，即可遍历2个像素点在二维离散坐标系下的“连线”所对应的像素点，获得连线的光强变化. 另外，之所以不以像素点的8邻域进行判断是由于使用8邻域判断出来的两点间“连线”的长度与实际长度也有较大的误差，这也是离散坐标系的对于长度的定义与连续坐标系的差异之处.

图7 步行进法示意图(蓝色线为已经走过的像素点)

4.2 线扫描光强变化拟合法

获取了加载点到圆心连线上的灰度变化如图8所示，连线上像素点到圆心的距离d=R-r，可以观察到，灰度变化曲线起伏成峰，在接近加载点位置(左)峰比较密集，在靠近颗粒圆心的位置，峰较为稀疏，实际上这些峰正是反映了连线上各级应力条纹的位置及条纹宽度. 这样，可以通过与理论对比的结果来反推加载力的大小. 对比实验中加载点到圆心连线灰度变化的峰的位置与理论的峰的位置，可以观察知道由于在峰谷与峰顶位置灰度变化不够平滑，原因是存在噪声的起伏. 由此，以较为平滑的腰的位置来确定峰的位置.图片为0～255的灰度级图片，若以I=Ia/2的光强来决定腰的位置，即以灰度级127确定腰的位置. 另外，灰度级变化以像素点位置离散，如一像素点灰度级为123，下一像素点灰度级为130，在这2个像素点之间，要找灰度级为127对应的位置则需以插值法进行插值找根，由于实验图像精度足够，2个像素点之间直接采用线性插值即可.

图8 某一加载点至圆心连线的灰度变化

值得注意的是，这样处理往往只能得到连线上靠近圆心的一些级数较低的峰的腰的位置. 这是因为实验中由于颗粒模型存在厚度，发生光的折射，加载点附近的峰的暗条纹不够明显，以及拍摄角度问题，在颗粒边界光强附近会有较大的误差. 但同样由于在接触点附近的条纹均为高级条纹，且理论的灰度变化在靠近加载点附近也奇异，对寻找接触力意义不大. 在处理时往往只处理低级的条纹，这是可以接受的，且无法避免. 如此获取了实验的加载点到圆心连线灰度变化的峰的腰的位置后，这些位置相当1组描述峰的腰的位置的数组向量. 有几个加载点就会有几个这样的向量，将其定义为s(i)，i表示对应的加载点. 同样的，由颗粒的多点加载光弹性理论，在加载情况下，也可从理论上得到描述加载点到圆心的连线灰度峰的腰的位置的向量. 理论得到的向量的维数往往比实验得到的向量的维数要多，这是由于由理论可直接得到更多高级的条纹峰的腰的位置. 但是在对比理论和实验的峰腰位置向量时，需要将理论的向量的维数截取到与实验的向量维数一致，相当舍弃掉了高级条纹，将这样处理后得到的峰腰位置向量定义为w(i,Fi)，i代表对应的加载点，Fi为第i个加载点的加载力.

要反推出加载点的加载力大小，接下来让加载点的加载力Fi在某范围内变化，计算出每种加载力大小的峰腰的位置数组向量，即计算出一系列的w(i,Fi)，减去实验对应加载点的峰腰位置向量w(i)，再取其模方，即得到了第i个加载点理论的加载力与实验的差值的平方. 而每个加载点在某一加载力下都能得到这样的差值方，对于一种加载情况，以加载点的加载力大小为索引，可建立对应不同加载情况理论与实验的差值方的i个维度矩阵C，其矩阵的矩阵元表示为C(F1,F2,F3,…，Fi)，i代表第i个加载点，Fi为第i个加载点的加载力，对其数学定义为

(10)

比如加载点数3，每个加载点的加载力的变化以1N为步长，且F∈[10,30]，那么矩阵C为20×20×20的3维矩阵，这个矩阵的最小元素值对应的索引(F1,F2,F3)则为需要的接触力的大小. 使用矩阵C寻找加载力大小的好处是它对于误差的包容性，因为理论每个加载力情况的结果与真实实验结果都存在相同的误差因子，误差因子包括：对低级条纹峰的腰位置的错误获取，颗粒应力光图的灰度级范围不在0～255. 这些误差因子有时难以甚至不能抹除，而矩阵C最小元素即对应误差因子存在时与实验的差别最小加载情况，也就是所需的最优结果.

通过多点加载二维颗粒实验，利用线扫描光强变化方法寻找加载力大小的过程如图9所示.

图9 线扫描方法流程图

如图10所示，以3个加载点为例，实验测得加载力与线扫描光强方法处理得到的加载力对比如表1所示.

(a)线扫描法获取连线上灰度变化的示意图

(b)线扫描法得到的结果模拟图10 线扫描法获得连线上灰度变化及模拟结果

表1 3点加载力对比表

将线扫描方法处理得到的加载力大小反过来模拟出颗粒的应力光图,如图10(b)所示，可以看出与实验图案图10(a)相比还是有所区别，这是由于实际实验中模型受加载后发生形变，加载处的加载方式已经由线接触变为了面接触，这也是通过先扫描的方式来寻找加载力的大小的原因. 如果对颗粒上全场的光强来考虑，一是由于3维圆柱颗粒的线与面接触方式不同带来的应力光图的差异带来的误差在考虑全场的光强时会更大；二是颗粒边界处的光强由于折射拍摄角度等带来的误差也难以处理；三是计算量较大. 但考虑颗粒上的全场光强的方法也有一定好处，即可以处理更加复杂加载情况下的图案. 要靠这种方式来寻找加载力大小，先让计算机进行预处理再人工调节寻找加载力大小，通过计算机与人的相互协作的方式反而更简单.

5 结束语

基于颗粒应力光图的灰度梯度平方分布图像的边界模板法可以获取颗粒的接触点位置. 该方法优点是能够在光弹颗粒条纹不十分密集的情况下，找出在颗粒边界处有较大的加载力存在的加载点的位置. 其缺点是不能够适用过于密集的颗粒应力条纹光图，且处理得到的加载点位置可能存在1°～4°的不确定度. 如果在大量颗粒接触时，可采用该方法与阈值切割方法[12]结合的思路. 在离散坐标系下寻找2点之间“连线”的算法，能够更加优化地处理Bresenham算法得到的“连线”的像素点个数在2点斜率在为非零有限值时不合适的问题. 基于对颗粒应力光图的线扫描的获取接触力大小方法，结合了二维颗粒多点接触的光弹性理论、数字图像处理，能够得到颗粒每一接触点的较为精确的接触力的大小. 为了提高该方法的运行速度与处理精度，应该事先确定颗粒上的加载力合理的模拟的范围，可采用平均灰度梯度方的方法得到颗粒的平均加载力以确定加载力的范围.