裂缝是沥青混凝土路面病害中最为常见的类型之一，采集路面影像数据可以真实地反应路面破损情况，它是路面养护中的一项重要环节。采用计算机进行自动裂纹检测可以更为客观合理地完成识别，从而避免仅依赖于主观传统检验方法带来的不足。

Guan等人^[1]通过综合运用阈值、增强及形态学腐蚀方法获得GRF(geo-referenced feature)特征，用于识别裂缝，并在文献^[2]中采用ITV(iterative tensor voting) 模型分类MLS(mobile laser scanning)数据，以此改善裂缝识别的准确性。文献^[3-4]采用小波变换及其变体，实现了裂纹提取技术。基于动态优化的方法^[5-7]可以有效处理模糊和不连续的路面图像，然而由于这种方法需要密集计算，因此降低了该方法的实用性。文献^[8]提出一种基于马尔可夫随机场的多尺度检测裂缝的方法，采用一维高斯平滑滤波器增强裂缝，并通过2维匹配滤波的方法完成检测。曹建农等人^[9]以分开—合并的方式通过Mean Shift算法完成图像分割，而后提取裂缝骨架，通过骨架内插值得到完整裂缝，实现了裂缝的识别与测量。文献^[10-12]采用机器学习及神经网络的方法完成了裂缝识别,但参数的选择依赖于裂纹的变化与图像质量，且需要大量样本数据进行训练。例如，文献^[12]以样本的学习为基础，通过子集图像数据的自动选择建立无监督训练系统，该系统可以将非重叠图像块分为包含与未包含裂缝的区域。

尽管上述方法可以较好地完成裂缝识别，但由于车辆在行驶过程中，周围环境将时刻发生变化，因此相机采集到的影像中可能存在多种干扰，具体包括：采集设备的成像质量，树木及车辆的投影，光照变化，油渍、树枝与稻草等条状物，各类垃圾等，而上述方法均未考虑到这类干扰对裂缝识别准确率的影响。例如采用文献^[8]的方法，会将树枝与稻草等条状物识别为裂缝，在这种情况下，裂纹的误识别率将会增加，从而降低路面病害估计的准确性，间接影响路面的养护工作。针对现有方法存在的问题，提出一种面向线阵列CCD(charge-coupled device)道路影像的裂缝鲁棒识别方法，主要贡献点为：

1) 通过分析干扰物与道路裂缝在线阵列CCD成像中呈现的特点，提出4种约束条件识别裂缝；

2) 分别对上述约束条件给出相应的数学求解算法，通过建立道路裂缝鲁棒识别流程确定裂缝位置，并给出裂缝的提取方法。

为了提高道路裂缝的识别率，对路面上可能影响裂缝识别的干扰物进行归类，图 1给出了各种类型的干扰物。图 1(a)为含有树枝与稻草的路面影像，与图 2中的裂缝相比，其形状特征相似，都可以归纳为线条形状；图 1(b)为路面油渍，除了部分区域形状与裂缝相似外，通常会伴有块状区域出现；图 1(c)为树木与车辆的投影，其投影区域与其他 区域的光照强度不同；图 1(d)中的环境更复杂些，包含多种类型的不同形状垃圾。

通过对比图 1与图 2中干扰物与裂缝的特点可见，这些干扰物均可能影响裂缝识别的准确率。尽管如此，通过线阵列CCD获取的道路影像仍可总结出以下特点：

1) 裂缝与沥青路面在局部区域中灰度值差异较大。文献^[12]对裂缝的形态特征进行了分析，其区域灰度值与邻域灰度存在较大差异，如图 2所示。

2) 裂缝在影像中呈现的连续性较差。该特点主要由线阵列相机采集精度不高的因素造成。

3) 裂缝的宽度与长度比值较小。与大部分的油渍、阴影等块状区域相比，裂缝的宽度远小于它的长度。

4) 同一段裂缝的走向基本一致。沥青路面裂缝类型通常分为横向裂缝、纵向裂缝、网状裂缝及块状裂缝，而后两类裂缝均由纵向与横向裂缝交错而成，因此本文将研究重点放在前两类裂缝的识别上，而横向与纵向的概念本身就具有一定的方向性，因此分割后的断裂区域也将在其方向上基本保持一致。

由于线阵列相机采用线扫描的方式，通常会产生光照不均匀的现象，因此采用分块的方式将待处理部分局部化，降低这一因素的干扰。合并所有分块后，获得对输入图像的初步分割结果。后续处理过程中，依据裂缝的4个特点，综合运用椭圆拟合方法及质心等区域特征，实现道路裂缝的识别，输出结果为裂缝所在的矩形区域与裂缝长度，详细的识别流程如图 3所示。

为了降低光照不均所产生的干扰，需要对图像进行分块，其分块大小需要依据采集图像的尺寸而定。通常由于采集设备是固定安置于车辆上的，其采集图像的角度及尺寸在整个行车过程中不变，因此在这种情况下，分块的大小只需设置一次。同时，块的尺寸不宜过大或过小，过大将增加光照的干扰，而过小时不符合本文提出的“裂缝在局部区域中占据较小的面积比”这一特点。例如，对于一幅4 096×2 048像素的道路影像来说，可以将分块尺寸设置为64×64像素。

从图 2的典型道路裂缝上看，裂缝区域与其他区域存在明显的灰度差异，因此选用水平集类方法中的CV模型对分块区域进行处理，原因在于该方法可自动地将图像分为灰度较高及灰度较低的两类区域。CV模型的“能量”泛函表达式为

$\begin{align} & E({{c}_{1}},{{c}_{2}},C)=\mu {{\int }_{C}}ds+{{\lambda }_{1}}\iint\limits_{{{\Omega }_{1}}}{{}}{{(I-{{c}_{1}})}^{2}}dxdy+ \\ & {{\lambda }_{2}}\iint\limits_{{{\Omega }_{2}}}{{}}{{(I-{{c}_{2}})}^{2}}dxdy \\ \end{align}$

(1)

式中，第1项为曲线C的全弧长，第2项与第3项分别是内部区域和外部区域的灰度值与标量c₁和c₂的平方误差。采用变分水平集法求解，在c₁和c₂固定的条件下，相对于u最小化能量泛函方程，通常取λ₁=λ₂=1，采用半隐式方案求解得到

$\begin{align} & u_{ij}^{n+1}=u_{ij}^{n}+\tau {{\delta }_{\varepsilon }}(u_{ij}^{n})[\mu Q(u_{ij}^{n+1})- \\ & {{({{I}_{ij}}-{{c}_{1}})}^{2}}+{{({{I}_{ij}}-{{c}_{2}})}^{2}}] \\ \end{align}$

(2)

$\begin{align} & Q({{u}^{n+1}})=D_{x}^{(-)}\left( \frac{D_{x}^{(+)}{{u}^{n+1}}}{\sqrt{{{(D_{x}^{(+)}{{u}^{n}})}^{2}}+{{(D_{y}^{(0)}{{u}^{n}})}^{2}}}} \right)+ \\ & D_{x}^{(-)}\left( \frac{D_{y}^{(+)}{{u}^{n+1}}}{\sqrt{{{({{D}^{(0)}}_{x}{{u}^{n}})}^{2}}+{{({{D}^{(+)}}_{y}{{u}^{n}})}^{2}}}} \right) \\ \end{align}$

式中，u为距离矩阵，τ为迭代步长系数，δ_ε(*)为正则化Heaviside函数的导数，D_*^(－)、D_*⁽⁰⁾ 、D_*⁽⁺⁾分别为前向差分、中心差分及后向差分。经典CV模型的收敛速度不稳定，初始轮廓的选取位置离目标越近，则收敛速度越快，反之越慢，因此降低了该方法的实用性。目前，已提出较多的改进方法^[13-14]解决这一问题。例如在文献^[13]中，对初始距离函数的构造进行改进，在欲分割的图像中，构造无穷多个圆形区域作为零水平级的初始区域，每个区域的半径无穷小，这种极限的初始化方法将使高1维度距离空间与零水平集平面重合，加速了CV模型的收敛速度。采用上述方法对道路影像进行二分类处理的目的是将局部块图像中的显著部分提取出来，获得二值图。

采用文献^[13]的方法对图 4(a)进行二分类，得到如图 4(b)所示的二值图。文献^{[9, 12]}对裂缝的形态特征进行了总结：裂缝区域的灰度值高于邻域灰度，因此将较高灰度的类置1，较低灰度类置0。采用区域增长算法标记每个置1区域，统计每种标记区域的面积，并根据经验阈值Tr_min及Tr_max过滤极小与极大面积区域，得到如图 4(c)的处理结果。

为了更好地完成过滤，采用形态学腐蚀算法以单位1为半径腐蚀一次，将可能存在“粘连”的小区域分开。尽管可能会将某些以单像素连接的裂缝截断，但由于本文方法是以判断分割区域共线为基础识别道路裂缝的，因此这种处理并不会影响其正确性。

道路裂缝具有一定的宽度与方向性，为了利用第1节给出的3、4项特征，拟采用结合椭圆拟合的方法^[15]求取每个标记区域的宽度与正方向。

如图 5所示，令S为待拟合区域，对S进行椭圆拟合，求解椭圆参数(包括长半径、短半径、长半径与水平方向的夹角、离心率)，建立以长半径为x轴、短半径为y轴、原点为o的坐标系。

根据椭圆离心率与阈值Tr_e去除非条状区域，并求取S的质心c。由于质心可能不在当前曲线上，因此以椭圆y轴为主方向，构造过c点的直线L，确定L与S的交点p，并统计p沿L方向的像素点数量，以此作为S的宽度值w，将w与给定阈值Tr_d进行比较，移除较大宽度的区域，实现对干扰的进一步过滤，结果如图 6所示。

依据道路裂缝的特点分析，结合图 2给出的几种方向的裂缝情况，将所有过滤后的区域依据椭圆x轴与水平方向的夹角θ进行分类，即

$Type=\left\{ \begin{align} & \begin{matrix} 1 & 0{}^\circ \le \theta <25{}^\circ 或155{}^\circ \le \theta \le 180{}^\circ \\ \end{matrix} \\ & \begin{matrix} 2 & 25{}^\circ \le \theta <75{}^\circ \\ \end{matrix} \\ & \begin{matrix} 3 & 75{}^\circ \le \theta <115{}^\circ \\ \end{matrix} \\ & \begin{matrix} 4 & 115{}^\circ \le \theta <155{}^\circ \\ \end{matrix} \\ \end{align} \right.$

(3)

分类的目的是将不同方向的线段归类，便于为下一步判断共线提供支持。式(3) 中，Type为分类号，不同的类别代表不同的方向。以水平方向分类为例，即Type为1时，根据影响裂缝鲁棒识别的第4项特点(同一条裂缝的方向基本一致)，对识别道路裂缝算法进行描述，过程如图 7所示。图 7(a)给出6条线段模拟裂缝，其中下半部分的4条线段属于裂缝部分，而上半部分作为干扰项。首先分别求取A-F每条线段的质心，并计算质心间的关系，用矢量表示如图 7(b)所示。

图 7(b)中表的基本单元描述了质心间的关系，用(矢量角，矢量模值)表示，如表中第2行、第3列的(355,10) 表示以质心A为原点建立坐标系，B在该坐标系中与水平轴的夹角为355°，与原点的距离为10个像素单位。从A与B的空间位置关系上考虑，A与B间的关系(355,10) 及B与A间的关系(170,10) 应是相同的，因此对图 7(b)表中180359间的值进行修正，将其映射到0179的空间上，得到图 7(c)的结果。令距离阈值Tr_dist=30，夹角范围取0°≤θ<25°或155°<θ≤180°，对图 7(c)中每个单元格中不符合该条件的值进行剔除，即表中红色标记的部分。最后统计表格中每行中列的数量n，该值将决定有多少条分割区域共线，共线数量越多，则判定为裂缝的概率越大。例如本例中令n≥2，即可将A-D与E-F分离，达到道路裂缝识别的目的。

上述方法存在的问题为：凭借当前行中1的数量判断是否为裂缝，在图 7(e)的情况下将出现误识别。尽管A与B、C、D的夹角符合0°≤θ<25°或155°<θ≤180°取值范围，但它们并不共线。为了解决该问题，给出递归方法处理矢量表中数据，其思想为继续判断列中值为1节点与其他节点间的关系，以此修正相关表项，伪代码如下：

Correct (&Table, &NewTable) //&表示对象引用

{m=size(Table); //通过size(*)求行、列数m

NewTable(1:m,1:m)=-1;

Tag=1; //裂缝编号

for i=1:m

{

  Mark (Table,NewTable,m,i,Tag);

  Tag++;

 }

}

bool Mark (&Table ,&NewTable,m,i,Tag)

{//sumones(*)用于统计1的个数

if (sumones (Table(i,:))<Tr_n) return 0;

for j=1:m

{

 if (NewTable (i,j) !=-1) return 1;

  NewTable (i,j) = 0;

 if (Table (i,j) = =1)

  if (Mark(Table,NewTable,j,Tag))

   NewTable (i,j)=Tag;

return 1;

}

}

Correct函数的参数有两个：Table为2维矢量表(如图 7(d)所示)，NewTable为修正后的新表。Mark函数包括4个参数，除Table与NewTable外，第35个参数分别为表格列数、行号及裂缝编号。Correct函数通过Mark递归识别裂缝，不仅可以解决图 7(d)中的问题，同时可以将不同裂缝线段标记出来。采用上述方法对图 6进行处理，得到如图 8(a)所示的标记结果，图 8(b)为根据标记结果在原图中绘制的最小矩形包围。

在道路裂缝损坏的统计中，可将图 8(b)中矩形对角线长度作为该图中的裂缝长度，通过这种方法可以快速地近似估计道路裂缝的长度，更适用于粗略统计长距离道路损害等级。

在需要对道路裂缝进行更为细致的分级评估时，可以将图 8(b)的矩形作为活动轮廓模型的初始嵌入函数进行真实边缘的求解，该方法的优点是曲线在进化过程中将保持封闭与连通。然而由于这类水平集方法的求解速度及结果受初始封闭轮廓影响较大，迭代过程中也可能会受到其他部分的干扰，因此给出以下方法提取裂缝：

1) 断裂分割区域连接。以图 9(a)为例进行说明，该图中包含3条线段S₁-S₃，采用第3.2节给出的方法得到p₁-p₃。确定S₁-S₃椭圆拟合x轴方向，并分别延x轴正负方向上计算距p₁-p₃最远两点作为线段端点。m_1.1与m_1.2分别为S₁的两端的节点，由于S₁与S₂间的真实边缘的方向未知，因此计算m_1.2与m_2.1的中心点m_1-2，并以m_1.2与m_2.1的节点间距离为直径画圆，以1填充其内部区域。同理，m_2-3的计算与其相同，经过这样的处理，得到如图 9(b)所示的结果。

2) 活动轮廓模型迭代求解。以断裂分割区域的连接结果为依据，确定区域的初始迭代边缘，将其作为活动轮廓模型嵌入函数的初值。采用文献^[16]给出的标准最小化PDE模型进行迭代求解，曲线在进化过程中将保持封闭与连通，即

$\frac{\partial u}{\partial t}=div(g\frac{\nabla ~u}{\left| \nabla ~u \right|})+\lambda \frac{u-f}{\left| u-f \right|}$

(4)

$g\left( r \right)=\frac{1}{1+\left( r/K \right)}$

(5)

$r=\frac{{{f}_{xx}}f_{y}^{2}-2{{f}_{x}}~{{f}_{y}}+{{f}_{yy}}~f_{x}^{2}}{{{(f_{x}^{2}+f_{y}^{2})}^{\frac{3}{2}}}}$

(6)

式中，u为嵌入函数，f为待处理图像，λ为大于0的比例系数，g为边缘函数，f_x、 f_y为一阶差分，f_xx、 f_yy为二阶差分。式(3) 给出一种g函数的具体形式，其中K为选定的常数，用于控制g的下降速率，如图 10所示。通过调节反差参数K，可以获得不同的曲率调节度，获得理想的阻力曲线。采用水平集的方法求解，对应的嵌入函数u的偏微分方程为

$\begin{align} & \partial {{u}^{n+1}}=\partial {{u}^{n}}+\Delta tg\left\{ (\sqrt{{{({{D}^{0}}_{x}~f)}^{2}}+{{({{D}^{0}}_{y}~f)}^{2}})} \right.\cdot \\ & (D_{x}^{-}N_{x,\varepsilon 1}^{+,n}+D_{y}^{-}N_{y,\varepsilon 1}^{+,n})+\text{ }{{(D_{x}^{0}g)}_{\Lambda }}\cdot N_{x,\varepsilon 1}^{-,n}+{{(D_{y}^{0}g)}_{V}}\cdot \\ & N_{x,\varepsilon 1}^{+,n}+{{(D_{y}^{0}g)}_{\Lambda }}\cdot N_{y,\varepsilon 1}^{-,n}+{{(D_{y}^{0}g)}_{V}}\cdot N_{y,\varepsilon 1}^{+,n}+ \\ & \left. \lambda \frac{{{u}^{n}}-f}{{{({{u}^{n}}-f)}^{2}}+{{\varepsilon }_{2}}} \right\} \\ \end{align}$

(7)

式中，D_x⁰ f、D_x^-f分别为水平轴上的中心和前向的差分，N_y,ε1^±,n为u水平轴上的估计值，(·)_Λ:=max(·,0) ，(·)_V:=min(·,0) ，Δt为迭代步长，ε₁与ε₂为极小正常数，∂uⁿ为当前时刻的嵌入函数，∂uⁿ⁺¹为下一时刻的嵌入函数。

根据方法1) 构造嵌入函数的初始值u_ij⁰，初始化式(4) 的u_ij⁰求解为

${{u}^{0}}_{ij}=\left\{ \begin{align} & \begin{matrix} 1 & {{f}_{ij}}==0 \\ \end{matrix} \\ & \begin{matrix} -1 & {{f}_{ij}}>0 \\ \end{matrix} \\ \end{align} \right.\left( i\in m,j\in n \right)$

(8)

由于这种方法已经使迭代边缘位于原图边缘附近，因此可以通过较少次迭代完成边缘提取。采用该方法处理图 4(a)，在原始图像中进行迭代求解，裂缝提取结果如图 11所示，最后将提取结果中的不同区域中心线长度作为裂缝长度输出。

采用CPU主频3.2 GHz、内存2 GB作为实验环境，MATLAB编程。采用1 000幅尺寸为4 096×2 048的道路影像进行实验，涉及的主要参数取值分别为Tr_min=10、Tr_max=300、Tr_e>0.97、Tr_dist=30、Tr_d=10、Tr_n=4，其中3.4节方法涉及的相关参数依据文献^[16]选取λ=0.1、K=10、Δt=0.1、迭代次数设定为50次。

选用图 1道路影像数据作为测试集，采用3.1节、3.2节方法对这些影像进行处理，得到如图 12所示的初始分割结果。

由实验结果可见，这部分区域中，大多数的形状特征不符合条状物特征，因此可以利用椭圆拟合结果剔除大部分干扰，结果如图 12第2行影像所示。在这部分影像中，与道路裂缝相比，油渍与垃圾的线条宽度较大，可以利用这种特点将其滤除。稻草与树枝检测结果与裂缝极其相似，都呈细条状，通过第1节给出的影响裂缝鲁棒识别分析中第2个特点进行区别，利用3.3节给出的方法将其剔除。

图 13给出了一组裂缝识别结果。图 13(a)为5幅包含横纹、纵纹的影像。采用本文给出的3.1节、3.2节方法对这些影像进行处理，得到图 13(b)。由该二值图可见，裂缝与沥青路面、油渍、阴影间的区别已较为清晰地呈现出来。利用本文3.3节方法对图 13(b)进行处理，得到如图 13(c)所示的结果。由图 13(c)可见，尽管存在部分倾斜幅度较大的裂缝区段，依然可以对这类进行较好地标记，验证了3.3节所给出方法的有效性。依据图 13(c)给出的标记线段，通过3.4节方法的处理获得如图 13(d)的最终裂缝提取结果。

从2 000幅道路影像中挑选包含道路裂缝的影像100幅，并按序号等间隔分别取出5组未含有裂缝的影像100幅，每组200幅组成数据集进行测试。由于检测结果中可能对含有裂缝及未含裂缝的影像产生误识别，采用分类指标统计的方法进行算法性能评测。

将本文方法看做一个二分类器，即将每组数据分为两类，计为正例和负例，分别为：

1) TP：被正确地划分为正例的个数，即实际为正例且被分类器划分为正例的实例数(样本数)；

2) FP：被错误地划分为正例的个数，即实际为负例但被分类器划分为正例的实例数；

3) FN：被错误地划分为负例的个数，即实际为正例但被分类器划分为负例的实例数；

4) TN：被正确地划分为负例的个数，即实际为负例且被分类器划分为负例的实例数。

首先计算表 1各类指标，以此为基础计算结果如表 2所示。其中正确率(Ac)表示被分对的样本数除以所有的样本数，该值越高则正确率越高，识别结果越好。灵敏度(Se)、特效度(Sp)及精度(Pr)分别代表对正例的识别能力、对负例的识别能力、精确性的度量。从各项指标的均值上看，均达到95%以上，验证了本文方法的有效性。图 14给出了测量与实际估计的长度对比，由图可见，采用矩形对角线作为道路裂缝的长度与实测长度存在偏差，而采用结合活动轮廓模型的提取结果更加接近裂缝实测数值。

算法的执行时间包括2个部分，其中检测时间不包括3.4节给出的裂缝提取时间。从实验结果上看，平均每幅影像的执行时间约1 min，而裂缝的提取时间约10 s，如采用多核CPU编程运行将缩减这一执行时间。

在实际应用中，结合数据采集车辆给出的地理坐标，通过这种方式估计需要养护路段的信息。由于现有方法均未考虑道路裂缝识别过程中的干扰问题，因此本文未给出与其他文献的实验结果的对比数据。

$\left\{ \begin{align} & Ac=\frac{TP+TN}{P+N} \\ & Se=\frac{TP}{P} \\ & Sp=\frac{TN}{N} \\ & Pr=\frac{TP}{TP+FP} \\ \end{align} \right.$

(9)

本文重点分析了干扰道路裂缝识别的多种因素，并提出采用4种约束进行道路裂缝检测，针对每种约束设计了相应的处理算法，重点给出了分割区域的分类、干扰过滤和利用分类结果识别裂缝的方法。实验分别从干扰物过滤、检测准确性及裂缝提取结果上进行了分析，弥补了现有方法在识别道路裂缝过程中所忽略的干扰物问题，实验结果验证了本文方法的有效性与实用性。从算法性能分析上看，本文方法仍存在道路裂缝误识别的问题，下一步将结合深度学习提高识别的正确率。