车辙变形是目前沥青路面主要的病害之一。矿质混合料的组成对其影响显著，因此陆续有研究者尝试通过多种力学试验，希望建立材料组成与力学性能的联系。但试验结果变异性很大，真正的规律淹没在噪声中，问题的根源在于材料内部颗粒群复杂的结构形态，只有深入到材料内部，即从细观角度出发，充分挖掘颗粒群结构的形态特征和变形特征，才能深刻揭示由材料组成引起的复杂结构形态对力学性能的影响机理，建立材料组成与力学性能的联系。材料的细观特征(接触形态^[1-3]、空间分布均匀性^[4-5]，集料长轴角度^[5]等)，及其与宏观变形^[6-9]的相关性或宏观变形预测方面^[10-13]的研究较多，但由于缺乏有效的颗粒群变形探测手段，在材料组成与颗粒群变形特征方面研究难以深入，因此，探索多目标、精度高、抗干扰能力强的混合料颗粒群变形跟踪方法十分必要。

目前在混合料颗粒群变形跟踪方面主要有3类做法。一是利用特征参数总体的变化描述内部结构的改变，如吴文亮^[14]通过计算车辙试件切片图像在加载前、后中颗粒角度分布的变化，得到荷载作用下的集料转动规律，该方法由于没有对集料个体在加载前后进行辨识，无法得到集料个体的变形或转动。二是利用图像相关法计算局部变形^[15-16]，但集料颗粒发生较大位移或转动时，计算精度无法保障。三是乔英娟^[17]捕获贴于试件表面的方块靶标角点，从而检测出颗粒群的变形，但该方法靶标特征过于简单，仅仅依靠方块的4个角点特征，而集料边界存在大量角点，靶标角点容易淹没于噪声中，从而导致自动化识别程度大大降低。本研究将方格阵列靶标印制在集料表面，通过跟踪靶标的位移、旋转实现颗粒上任意一点的变形计算，而且组成靶标的角点阵列具有明显的几何特征，便于计算机进行识别。同时，角点阵列提供多个亚像素坐标用于变形计算，确保了精度问题，为解决混合料内部颗粒群变形跟踪难题提供了有效途径。

方格靶标由黑白相接的方格组成，两个黑色方格相交处具有明显的角点特征，称为方格靶标角点，如图 1(a)中的方格靶标形成4×4角点阵列。当靶标制成印章，印制在集料表面上时，对于变形前图像从上到下逐行扫描到的角点阵列中，第一个角点定义为角点1，其他角点序列的编号如图 1(b)所示。

图 1 方格角点阵列靶标 Fig. 1 Grid corner array target

基于方格靶标的颗粒群变形跟踪方法就是利用方格靶标的色彩、几何特征从变形前、后图像中的集料纹理、边界、划擦、色彩蜕变等干扰噪声中提取出来，再根据靶标角点阵列的几何特征实现集料颗粒位移、旋转等变形参数的计算。

方格靶标角点阵列的几何中心定义为靶标中心(x_t, y_t)，其变形后坐标为(x′_t, y′_t)。过角点1和靶标中心的直线为靶标u轴，过角点13和靶标中线的直线为v轴，两轴的方向如图 2所示。定义靶标v轴正半轴离开x正半轴逆时针旋转的角度为靶标的角度。将方格靶标制成印章，使用光敏印油将其印制在集料颗粒表面，颗粒发生位移或转动后，通过跟踪靶标中心的位移和靶标角度的转角可以计算出集料上任意一点的位移量，如图 2所示，靶标初始角度为a，变形后角度为a′，集料上任意点为(x, y)，则变形后的坐标(x′, y′)为：

(1)

图 2 变形跟踪原理 Fig. 2 Deformation tracking principle

使用方格靶标对混合料颗粒群进行变形跟踪，首先将混合料进行切割、打磨，使截面集料颗粒表面尽量光滑、平整，然后将靶标印制在集料颗粒表面上，风干后采集颗粒群初试位置图像，然后将试件组装、加载，采集加载后的颗粒群图像，通过捕捉两张图像中的靶标，获得各靶标的中心点及靶标角度，即可按式(1)计算集料上各点的位移情况，制作颗粒群的变形图谱。因此，方格靶标识别成为颗粒群变形跟踪的关键问题，其中主要难点为靶标角点阵列的定位和靶标几何参数计算。

图 3为带有靶标的颗粒边界的提取情形。

图 3 带有靶标的颗粒边界提取 Fig. 3 Extracting particle boundary printed with target

印制在颗粒群表面上的靶标识别为多目标识别，它是以单靶标逐个识别的方式实现的。为排除颗粒边界的影响，使用多阶段自适应边界分割方法^[18]对试件截面图像进行二值化。接着对各颗粒进行边界跟踪，对跟踪到的颗粒进行编号，要求每个颗粒有唯一编号。然后使用该编号对颗粒区域像素进行赋值，制成颗粒区域的编号映射图，并对映射图做腐蚀运算。由此，在利用外接矩形进行颗粒单靶标图像识别时，取映射图中与当前颗粒编号相同的像素区域，这样就将其他颗粒及颗粒自身的边缘影响排除了。

当靶标角度为kπ/2+π/4时，靶标方格的4边分别与图像的水平、垂直方向平行，称为平行靶标，它是靶标位姿的一种特殊形态，此时靶标图像具有明显的几何特征。首先使用Sobel边缘检测法对靶标图像进行边缘检测，排除颗粒边界的边缘像素，得到图 4(a)所示的边缘强度图像，并对其进行二值化，使边缘位置像素为0，非边缘区域像素为255。

图 4 靶标角度为kπ/2+π/4时的图像特征 Fig. 4 Image feature when target angle is kπ/2+π/4

对边缘强度图像做水平、垂直方向的投影：

(2)

(3)

式中，T_h(i)为第i行水平方向投影；T_v(j)为第j列垂直方向投影；m，n为颗粒图像长和宽尺寸；E(i, j)为边缘图像中第i行、第j列的边缘强度。

如图 4(a)所示，水平、垂直两个方向的投影均出现若干峰值，峰值的位置对应方格的竖向边界V₁~V₄和水平边界H₁~H₄，部分噪声在两个方向的投影曲线中也出现若干峰值，但其大小一般小于方格边界处的投影峰值。为了找到方格边界所对应的峰值位置，以方格边长r为周期做两投影的累加曲线C_h和C_v：

(4)

式中，λ为图像水平投影位置；k为方格从上到下的递增序数。

(5)

对于C_h，k为0，1, …, m/r；对于C_v，k为0，1, …, n/r。C_h，C_v曲线峰值位置分别对应着靶标方格的4条横、纵边界线。设ξ和ζ分别使曲线C_h，C_v达到峰值，通过搜索曲线C_h和C_v的峰值位置，找到ξ和ζ，方格靶标纵向4条边界V₁~V₄对应T_v曲线kr+ξ中4个最强的峰值位置，而方格靶标横向4条边界H₁~H₄对应T_h曲线kr+ζ中4个最强的峰值，也即由累加曲线的ξ和ζ可以重构平行靶标的16个角点坐标，如图 4(b)所示。

由上述分析，当颗粒靶标角度为kπ/2+π/4时，可通过靶标明显的纵、横投影曲线及各自以r为周期的累加曲线峰值位置计算平行靶标角点阵列坐标。

为找到任意位姿的靶标角点阵列坐标，对颗粒靶标图像逆时针旋转1周，得到每个旋转角度下的图像边缘强度图像，由边缘强度图像计算边缘强度的纵、横投影及相应的以r为周期的累加曲线。搜索各旋转角度θ下的边缘强度累加曲线峰值C_h(ξ_θ，θ)和C_v(ζ_θ，θ)，其中ξ_θ为靶标旋转角度为θ时，使纵向边缘强度投影以r为周期的累加曲线C_h取得峰值的点；ζ_θ为靶标旋转角度为θ时，使纵向边缘强度投影以r为周期的累加曲线C_v取得峰值的点，绘制C_h，C_v与θ的关系图，见图 5。

图 5 旋转角度与边缘强度投影累加峰值的关系 Fig. 5 Relationship between rotation angle and edge strength projection accumulation peak

从图 5中可以看出，曲线C_h(ξ_θ，θ)与C_v(ζ_θ，θ)的形状相同，相位差为π/2，靶标旋转1周的过程中，两曲线均出现4个明显的强峰值，且峰值位置相同。很明显，这4个峰值位置正是靶标角度为kπ/2+π/4时，方格的边界线分别与图像的水平、垂直方向平行的情况。由于角点阵列是关于靶标中心的中心对称，在计算靶标角点阵列坐标时，可通过搜索曲线C_h(ξ_θ，θ)或C_v(ζ_θ，θ)在[0, π/2]上的最大峰值计算靶标角度。设图像逆时针旋转角度为θ_m时，曲线C_h(ξ_θ，θ)在[0, π/2]为最强峰值，则靶标角度必然为θ_m，π/2+θ_m，π+θ_m，3π/2+θ_m中的一个。图像在逆时针旋转θ_m后，得到ξ_θ和ζ_θ值。在T_h曲线的kr+ξ上取4个最强峰的位置，即为V₁~V₄的横坐标；在T_v曲线上的kr+ζ上取4个最强峰的位置，即为H₁~H₄的纵坐标位置，由此可计算得到16个靶标角点在图像旋转θ_m后的坐标。最后，将这些坐标值顺时针旋回θ_m后得到靶标角点坐标，如图 6所示。

图 6 获取靶标角点阵列坐标 Fig. 6 Obtaining target corner array coordinates

得到靶标角点阵列坐标后，为了进一步提高坐标计算精度，以这些坐标为初始位置，对各角点进行Harris角点亚像素搜索。该算法根据角点在任意方向移动均会发生较大灰度变化这一特点，通过在高斯平滑后的图像上移动1个小窗口，根据小窗图像差分及实对称性质，构建参数R_corner用以表征窗口内的灰度变化与由图像梯度组成的二阶矩阵Q相关性。R_corner越大则当前位置是角点的可能性越大，R_corner由下式计算：

(6)

式中，；I_x和I_y分别为沿x和y方向的灰度梯度；det为矩阵的行列式；trace为矩阵的迹。

当R_corner在移动窗口(本研究取3×3窗口，步长为0.1像素)范围内最大时，即为角点，从而实现在给定角点初始位置的情况下的亚像素搜索。最后由16个亚像素角点坐标的均值获取靶标中心点坐标。

由前文分析可知，靶标角度可能为θ_m，π/2+θ_m，π+θ_m，3π/2+θ_m中的一个，但究竟是多大角度，仍需做进一步判断。对于变形前图像，根据角点编号的规则，角点阵列中的第一角点是图像从上到下逐行扫描过程中第一个检测到的角点，以该角点为中心取周边M×M的子图像S₀(M为奇数)，然后对原颗粒图像分别逆时针旋转θ_m，π/2+θ_m，π+θ_m，3π/2+θ_m，取各图像中以左上角角点为中心、大小为M×M的子图像S₁~S₄，计算S₀与S₁~S₄的相似度，其中相似度最大的图像所对应的旋转角度即为靶标角度，相似度计算式为：

(7)

式中，k为1~4，M取值为5。

对于变形后图像，设图像逆时针旋转角度为θ′_m时，曲线C_h(ξ_θ′，θ′)在[0, π/2]为最强峰值，对变形后的图像分别逆时针旋转θ′_m，π/2+θ′_m，π+θ′_m，3π/2+θ′_m后，取各图像中以左上角角点为中心、大小为M×M的子图像Q₁~Q₄，用式(7)计算S₀与Q₁~Q₄的相似度R_angle，该值最大的子图像所对应的旋转角度即为变形后图像靶标角度。

采用方格靶标跟踪颗粒群变形时，靶标的最小边长可以做到8 mm，这样的尺寸可以印制在9.5 mm粒级以上的颗粒表面上。尽管少量靶标图案超出集料颗粒边界仍可识别，但为了保证识别率，在印制靶标图案时，尽量将靶标印制在颗粒边界内部，且离边界越远越好。如图 7所示，将车辙试件切割后印上靶标，使用高分辨率扫描仪(300 dpi)采集变形前截面图像，然后将试件重新装入试模中进行车辙试验，加载完成后，取出试件，采集变形后的截面图像，分别对图 7(a)~(b)中的靶标进行捕捉，得到各靶标的中心点和靶标角度，按式(1)计算集料上栅格点的位移，制作变形图谱，见图 7 (c)。为了突出显示颗粒变形情况，将位移放大了3倍。

图 7 基于十字靶标的颗粒群变形跟踪 Fig. 7 Tracking particles deformation based on cross target

从图 7(a)~(b)中可以看出，在集料纹理、边界、表面划擦、色彩蜕变等干扰噪声下，几乎所有的方格靶标均被捕捉，尤其是车轮加载区域，颗粒变形大，存在较大的旋转角度，沥青从集料缝隙中挤出，对靶标图像造成了较严重的污染，在此情况下，方格靶标仍能被准确识别，说明方格靶标明显的几何特征及较为复杂的靶标图案具有较强的抗噪声干扰的能力，适用于大变形跟踪。

从图 7 (c)中可以看出，荷载附近区域的颗粒以向下运动为主，并沿着深度方向向两侧扩散。当这种扩散被左右边界试模约束时，颗粒位移的形态逐步变为向上隆起，但左右两侧的变形情况并不对称。在外力作用下，混合料内部以颗粒为单位不断改变运动方向，颗粒运动趋势受到混合料内部颗粒的空间排布、集料接触形态的影响。基于方格靶标的颗粒群变形跟踪作为一种有效的分析工具，将混合料的宏观变形与内部颗粒群的细观变形结合起来，并与反映混合料内部细观形态的集料形状、颗粒空间分布、颗粒接触状况等参数存在天然密切联系，因而颗粒群变形能够将宏观力学性能与细观结构形态参数联系起来，为深刻揭示混合料宏观行为的内在机理和改善混合料性能提供一种新的思路。

(1) 基于方格靶标的沥青混合料截面内部颗粒群变形跟踪方法，将靶标图案印制在混合料截面集料颗粒上，通过捕捉变形前、后图像中的靶标，获得其旋转角度和中心点坐标，实现集料颗粒上任意一点的变形跟踪。靶标图案边缘特征明显，结构较复杂，且提供了多个角点用于靶标中心的计算，具有抗干扰能力好、适用性强的优点，为探索沥青混合料内部颗粒群变形机理提供了一种有效的分析手段。

(2) 平行靶标图像具有明显的几何特征，它的边缘强度在水平、垂直方向的投影在方格边界处会出现强峰，结合其累加曲线峰值，可得到平行靶标角点阵列。根据靶标角度将平行靶标旋回可得到靶标角点阵列，最后通过靶标角点阵列亚像素均值化处理完成靶标中心计算。

(3) 由颗粒图像旋转角度与边缘强度投影累加峰值的关系可知，此曲线在[0, 2π]上会出现4个峰值，峰值相差π/2，且竖直方向与水平方向曲线形状相同，相位差为π/2，将4个峰位作为靶标角度的候选值，并结合变形前第一角点与各候选值下的平行靶标左上角角点的相似度分析，能够得到靶标角度。