吴 琼，薛 宁，刘春秀

（1.中国科学院空天信息创新研究院传感技术国家重点实验室（北方基地），北京 100190；2.中国科学院大学，北京 100049）

在冬季气温较低时，高速公路和机场路面结冰、积雪成为诱发交通事故、影响飞机安全起降的直接原因，精确监测路面冰雪情况并及时清除可以有效避免事故的发生。目前，国内外针对路面状态检测的研究主要以非接触式为主，如几何光学测量方法、近红外多光谱法[1]、光学偏振法[2]、雷达检测[3]和光强法等，但分别存在实际安装困难、检测点少、不能定量检测、成本高等不足。利用分类器的检测方法时，准确率受图像质量和分类器选择的影响[4]。

基于以上研究现状，根据在可见光范围内冰的高透射率和雪的高反射率[5]导致成像时灰度值不同，文中采用机器视觉方法对路面冰雪断面进行检测，断面通过已有模拟路面平台获得。与前述方法相比，该方法成本低且安装简单，可以定量检测冰雪厚度且无需大量气象数据支持。其难点在于自然结冰或积雪过程中冰雪厚度、边缘形状、结冰或积雪还是冰雪同时存在都具有随机性，传统的边缘检测算法已无法适用。文中首先通过机器视觉系统设计与测试保证量程和精度满足要求；其次，通过比较不同滤波器的去噪效果，采用两种保边滤波器和中值滤波器相结合的方式去除噪声；最后，设计了两种基于Canny 算子的有效边缘识别方法，实验测得两种方法的检测平均误差小于10%。

1 系统设计与算法实现

1.1 系统设计与分辨率测试

机器视觉系统性能与各部件均有关，需要根据实际应用需求进行系统设计[6]。根据现有产品，选用芯片尺寸为1/2.5 英寸，像素大小为2.2 μm×2.2 μm。根据奈奎斯特采样定理，设定芯片在像面处最小可分辨尺寸为8.8 μm，当视场范围和空间频率分别为40 mm×40 mm 和5 lp/mm 时，根据薄透镜放大率计算公式，如式（1）所示：

式中，h′和h分别为像高和物高，β为放大率，可以得到放大率应在0.088～0.107 5 之间。小孔成像中的比例关系如式（2）所示：

式中，L、h、f′、h′分别为物距、物高、焦距以及芯片高度。当物距为7～10 cm 时，得到镜头焦距在6.16～10.75 mm，文中选用8 mm 焦距。实验测得当物距为8.4 cm 时，放大率为0.088。最佳照明策略是最大化场景对比度，将重要区域的亮度置于相机的动态范围内，依据冰、雪在可见光范围内的高透过率、高反射率以及不同光源的应用特点[7]，文中采用白色环形LED 照明。

为了测试系统的实际空间分辨率，采用对比度法[8]在物距为8.4 cm 处利用不同空间频率的矩形靶标对不同视场下的调制传递函数MTF（Modulation Transfer Function）进行测试。MTF 计算式如式（3）所示：

式中，M′和M分别为输出和输入调制度，Imax、Imin和I′max、I′min分别是给定空间频率V在输入图像和输出图像上像素灰度的最大值和最小值，输入调制度M近似为1，含有系数π/4 是由于该处使用了矩形标靶而不是正弦标靶。利用程序求解I′max和I′min，通过计算得到系统MTF曲线如图1所示。以MTF大于0.6作为清晰成像的判断依据，比较图中曲线可以得到系统分辨率约为4 lp/mm，小于相机的理论值5 lp/mm，这是由于标靶精度、平整度及光源均匀性导致的。通过以上系统设计及测试，得到系统的实际测试范围和精度分别为0～40 mm 和0.125 mm。

图1 系统MTF曲线

1.2 冰雪厚度检测算法的实现

算法的基本流程如图2 所示，对读取到的图像进行对比度提升以便于后续识别，相机畸变校正中采用操作简单且精度较高的张氏标定法[9]，忽略高阶畸变系数和微小的薄棱镜畸变，选取合适的棋盘格[10]标定，得到相机内参和畸变系数如式（4）所示，重投影误差为0.38，将得到的相机参数值作为输入参数进行畸变校正。通过以上预处理后对图像进行去噪、边缘检测和提取。

图2 算法流程

1.2.1 图像去噪及Canny算子边缘检测

好的去噪效果是边缘检测的重要基础。冰雪图像噪声主要源于透光孔上的灰尘、水滴、冰中的气泡以及眩光效应，水滴和气泡可以利用统计排序滤波来减小。经过多次效果对比，文中采用非局部均值去噪、中值滤波和双边滤波结合的方式来去除噪声。图3 为非局部均值去噪的执行过程，以y为中心的邻域窗口在以x为中心的搜索窗口中滑动，计算两邻域窗口V(x)、V(y)之间的高斯加权欧氏距离来确定其相似程度，具体计算如式（5）～（7）所示：

图3 非局部均值去噪的执行过程

其中，(x)为经过滤波的图像，w(x,y)为邻域窗口y的权重系数，n为平滑参数，Z(x)为归一化系数。h的大小需要根据图像中的噪声水平来确定，在程序中该值设定为50。然后利用模板大小为27 的中值滤波器去除表面较小水滴引起的噪声点[11]，最后采用双边滤波器再次去噪。双边滤波器中的模板权值与像素间的欧氏距离和灰度相似性有关[12]，灰度差值和距离越大时权值越小。其空域模板系数、值域模板系数如式（8）和式（9）所示：

式中，(k,l)为模板窗口的中心坐标，(i,j)为窗口的其他像素坐标，σd、σr为高斯函数标准差，分别设定为15 和80，f(i,j)为待处理图像中坐标(i,j)处的像素灰度。

经过3 种滤波器对图4（a）降噪后得到的Canny算子边缘检测结果如图4（b）所示，检测到的边缘清晰且消除了由于气泡存在导致的冰块分层及灰尘和部分水滴的干扰。而当只使用其中两种滤波器去噪时，噪声较大，如图5（a）、（b）所示；图5（c）中的噪声也多于图4（b），这表明非局部均值去噪和中值滤波去除了大部分噪声，且未使边缘模糊而影响检测，3种滤波器结合的效果更好。

图4 经过3种滤波器去噪后的Canny算子检测结果

图5 不同去噪情况下的Canny算子检测结果

为了比较不同边缘检测算子[13]对冰雪图像的检测效果，依次采用Roberts 算子、Prewitt 算子、Sobel 算子和Laplacian算子对降噪之后的图4（a）图进行检测，结果如图6所示，从图中看到Roberts算子检测效果最差，Laplacian 算子未检测到冰雪分界，而Prewitt 算子和Sobel算子效果相似，虽然检测到了冰雪分界但比较模糊，为多像素宽度。图4（b）所示的Canny算子检测结果更为清晰，所以文中采用经典Canny算子，其基于信噪比、定位精度和单边缘响应准则，实现过程可分为4步。

图6 各算子边缘检测结果

1）对图像进行高斯平滑以提高信噪比；

2）采用如式（10）所示Sobel 算子与图像进行卷积得到梯度大小和方向，g(x,y)、α(x,y) 如式（11）、式（12）所示：

其中，Gx和Gy分别为水平和竖直方向的一阶导数值；

3）在梯度图内沿着每点的梯度方向判断该点邻域内的极大值，若该点不是极大值则将其置零，即进行非极大值抑制，保留局部梯度最大的点，以得到细化的边缘；

4）设定两个阈值T1、T2且T1＜T2，小于T1的点置零，大于T2的点赋最大值，在两者之间的像素点使用8 连通域确定其是否为边缘点，只有与大于T2的像素相连接时才被认为是边缘点[14]。与其他边缘检测算子相比，Canny 算子不仅得到了细化的边缘且更加准确。

1.2.2 两种有效边缘识别方法

经过以上3 种滤波器的去噪和Canny 算子进行边缘检测，可以看到所检测到的冰雪边缘不规则分布，且不同于一般工业零部件的厚度检测，冰雪厚度的可能值在0～20 mm 范围内变化，待测范围小且厚度值不固定，有效边缘定位成为整个厚度检测算法中的难点。为此，文中设计了基于边缘特征的一阶差分法和自适应阈值法，以实现有效的边缘检测。理想的水平边缘经过Canny算子检测，得到的边缘纵坐标曲线如图7（a）所示，但实际上边缘倾斜、不平整且有噪声影响，导致Canny 检测结果如图7（b）所示，可以看到数据量过大且曲线趋势复杂，为了在大量数据中精确定位边缘值，首先设计了基于边缘特征的一阶差分法。

图7 水平边缘各点纵坐标曲线

基于边缘特征的一阶差分法基本思想：利用相邻边缘坐标差值和有效边缘长度这两个边缘特征，得到边缘坐标数组中各有效边缘之间的分界位置，然后求解各段坐标均值，作为最终提取结果。其实现过程包括如下步骤。1）获得各点纵坐标数组：对Canny 算子检测结果进行概率霍夫变换[15]，得到被检测为边缘的各点纵坐标组成的数组，对数组中元素排序得到数组Y；2）滤除平滑边缘：利用式（13）计算数组Y的一阶差分ΔY，设定a=2，并与ΔY[i]进行比较，得到ΔY中大于2 的元素相应的索引数组A，有效边缘点的索引就包含在内；3）得到长度大于b的边缘：继续对A求一阶差分得到ΔA，设定b=50 并与ΔA作比较，得到ΔA中大于50 的元素相对应的索引数组C；4）求Y中阶跃变化的位置：根据各索引数组的对应关系，以C中索引为界，分段求解ΔY中最大值作为最终分段，识别Y中坐标值的分界数组D，求解各段坐标均值，即得到边缘提取结果。第一步中采用概率霍夫变换，是利用其统计特点来减小数据量。基于边缘特征的一阶差分法将直接求解分布范围较窄的阶跃点转化为先求解分布较宽、长度大于b的边缘位置，再据此求得阶跃位置。该方法不仅有效排除了噪声的干扰且检测精度和可检测到的边缘长度可以根据应用需要进行调节，其检测精度由相邻边缘点纵坐标差值a决定，可检测到的边缘长度由b决定。

自适应阈值法不同于上述基于边缘特征的一阶差分法，其是根据检测范围和精度来选取阈值数组，将阈值数组中每一个元素依次与坐标数组一阶差分的元素进行比较，从而得到适应于待测厚度的阈值和厚度个数，该阈值是消除噪声和干扰边缘对最终结果影响的关键，厚度个数作为冰雪识别的辅助条件。该方法可分为3 个步骤。

1）获得各边缘坐标Y：由霍夫变换[16]得到边缘线后，计算各线与图像左、右边界线的交点，两交点纵坐标p1、p2 计算式如式（14）所示：

其中，ρ为极坐标中直线与原点的距离，θ为直线的垂线与极轴的夹角，w为图像宽度，ρ和w的单位为像素；

2）求解符合当下检测厚度的阈值和厚度个数：将各边缘线纵坐标数组Y的一阶差分ΔY中各元素，依次与阈值数组Thresh[40,30,20,10,5,2]中的元素比较，得到的阈值作为识别Y中有效边缘坐标的依据，峰值个数表示在该阈值下检测到的厚度值个数；

3）识别待测物质：结合峰值个数和冰、雪灰度分界值150 对冰和雪进行识别，并通过相机放大率得到最终厚度。该方法的检测精度由Thresh数组中的最小值决定，根据不同检测范围和精度可具体设定数组中各元素值的大小。

ΔY和数组Thresh比较的过程是自适应阈值法的关键，主要可分为两步：第一步求解3 个索引数组Index1、Index2 和Index3，统计ΔY中大于Thresh[i]的元素个数，将其依次存于数组number0 中，并将number0 中元素值为1、2 和大于2 的相应索引分别存于数组Index1、Index2 和Index3 中，该索引值即为对应中阈值元素的索引。由于待测厚度和厚度值个数的随机性，所以第二步根据厚度检测中可能出现的情况将以上数组分为3 种情况进行判断。3 种情况包括：1）Index1 不为空；2）Index1 和Index2 都为空，此时表示超量程或厚度小于检测精度；3）Index1 为空，Index2 不为空，此时厚度值个数为2。Index1 不为空时，根据边缘数据特点，又分为Index2 不为空、Index2 和Index3 都为空、Index2 为空且Index3 不为空3 种情况。这种选取合适阈值的方法很好地解决了在待测物厚度、边缘形状和待测物个数随机变化的情况下对有效边缘的误判问题。

2 实验结果及分析

2.1 两种边缘检测算法验证

在环境光为室内灯光，物距为8.4 cm 时，采用聚丙烯酸钠制造的人造雪代替真实雪，在LED 照明下得到如图8（a）所示的冰雪图；去噪后利用Canny 算子检测边缘得到图8（b）所示的冰雪边缘，冰表面存在较多气泡引起的噪声，黑色底座因未在聚焦范围，所以边缘未被检测到；图8（c）为图8（b）中所有灰度值不为0 的像素点的纵坐标曲线，可以看到曲线没有明显的水平边缘规律性，且数据量很大。而经过两种方法处理得到的图9 中，数据量明显减少且曲线上圆点标记值即为利用基于边缘特征的一阶差分法和自适应阈值法得到的有效边缘均值计算结果，可以看到，结果都位于每段近似水平曲线的中间位置，表明了边缘提取结果的准确性。

图8 Canny算子检测结果

图9 两种有效边缘识别方法得到的结果

为了定量描述检测结果的准确性，如表1 所示为边缘检测结果、冰雪计算厚度、实际厚度和误差，误差计算如式（15）所示，Treal和Tcal分别为实际厚度和计算厚度值。实际厚度的测量以畸变校正后图像中央位置处对应的冰雪厚度为准。从表1 中可以看到两种方法的误差都小于10%。

表1 两种方法检测结果

为了避免偶然性，改变光照条件、边缘形状和冰雪厚度，对图10 的6 幅图片分别进行厚度检测，图10（a）、（b）、（c）分别为室内环境光下的成像，图10（d）、（e）为LED 灯照明下的成像，图10（f）利用LED 照明并且增加了实际应用时的透光板，引入了灰尘和眩光效应带来的干扰。相应检测结果如表2所示，可以看到两种方法对冰雪检测的平均误差都不超过10%。

图10 不同边缘形状和厚度的冰雪图像

2.2 误差分析

通过以上实验可以发现，基于边缘特征的一阶差分法整体误差高于自适应阈值法，这是由于该方法基本原理是基于水平边缘坐标特点的，当冰雪边缘不平整程度远大于检测精度时，会在一条边缘处出现多个阶跃位置，而算法只提取一阶差分最大位置。此外，在实验操作过程中，物距测量存在误差，导致实际放大率与预设的0.088 放大率不一致，从而在计算冰雪厚度时引入误差。

3 结论

文中提出使用机器视觉的方法对路面冰雪情况进行探测，相比于其他已有方法，其优点在于成本低，且可以定量检测冰雪情况。针对该方法在具体实现中的难点，文中主要进行了以下工作：

1）通过系统性能设计确保量程和精度满足检测需求；

2）针对冰雪图像噪声特点，采用两种保边滤波器和中值滤波结合的方式进行降噪处理；

3）针对自然结冰积雪过程中厚度范围、待测物个数和边缘形状的随机性，在Canny 算子的基础上，设计了基于边缘特征的一阶差分法和自适应阈值法提取边缘，两种方法的检测平均误差在10%以内。

以上工作为以后机器视觉应用于实际路面的冰雪厚度监测提供了可能。