一种隧洞衬砌混凝土管片表面粗糙度的测量方法与流程

1.本发明涉及混凝土材料表面图像处理应用领域，特别是一种隧洞衬砌混凝土管片表面粗糙度的测量方法。


背景技术：

2.隧洞衬砌混凝土管片被广泛应用泥水盾构隧洞工程中，是一种支持、维护隧洞长期稳定的永久结构物。其混凝土管片表面粗糙度是计算隧洞壁面糙率，衡量管片与周边岩土或二次衬砌混凝土胶结面强度与密封程度的重要参数。研究高精度、高效率的混凝土管片粗糙度测量方法具有重要的工程价值和理论意义。
3.目前，混凝土表面粗糙度的测量手段有很多，主要分接触式和非接触式两种类型。其中，接触式测量方法包括触针法、铺砂法、比较法、印模法等，该类方法是对表面粗糙程度的一个直接测量，操作简单直观。然而，这种测量方法一方面会对混凝土管片表面造成一定程度的损坏；此外还操作耗时，费力，不适于大面积进行测量。非接触式测量方法，大都涉及到光学仪器，如光切法，激光扫描法、全息法等，此类测量方法精度高，速度快，对待测面不产生损伤，但其测量设备购置成本高，操作复杂，当前主要用于实验室研究或大型、重要工程的测量，在中小型工程普及应用存在一定困难。
4.随着计算机的迅速发展，计算机视觉和图像处理技术得到快速普及，数字图像技术逐渐被工程测量领域接纳采用。该技术结合于图像测距原理，借助摄影设备便能完成图像采集，操作简单、成本低廉。已有研究中，数字图像技术已经被引入混凝土表面粗糙度的测量领域，但在实际工程中，数字图像技术还存在着如下不足，有待进行改进：
5.1、数字图像技术易受测量环境、操作手法的限制，无法获取优质、统一的图像数据，从而无法保证计算所得粗糙度的精度。
6.2、图像数据的数字化处理方法尚未实现简洁、易操作性，高效率、高精度的图像特征参数计算方法还有待研究。
7.2、现有的粗糙度计算公式多达上百种，引用的参数包括粗糙度轮廓指数、构造深度、起伏均方根等十余种，涉及露石混凝土表面，新、老混凝土表面等多种表面特点，且均取得了较好的应用，由此说明不同的结构表面存在其独特性，粗糙度的表征计算公式也存在差异，目前尚未出现针对衬砌表面粗糙度的成熟的计算公式，因此对其进行表征计算时，因站在已有的常见通式的基础下，结合真实条件，拟合得到高精度的计算公式。


技术实现要素：

8.本发明要解决的技术问题是针对上述现有技术的不足，而提供一种隧洞衬砌混凝土管片表面粗糙度的测量方法，该隧洞衬砌混凝土管片表面粗糙度的测量方法根据混凝土管片本身特点及测量条件，进行混凝土管片表面图像采集和处理，提取像素点灰度值等图像数据，引入图像特征参数值来描述混凝土管片图像粗糙度，因而能够得到高精度的混凝土管片表面粗糙度。
9.为解决上述技术问题，本发明采用的技术方案是：
10.一种隧洞衬砌混凝土管片表面粗糙度的测量方法，包括如下步骤。
11.步骤1、标定测量子面：测量子面包括计算子面和特征子面；具体标定方法包括如下步骤。
12.步骤1a、确定计算子面的大小和数量：根据待测量的隧洞衬砌混凝土管片表面整体面积，确定计算子面的大小和数量；计算子面为在待测量隧洞衬砌混凝土管片表面均匀分布的子面。
13.步骤1b、标定计算子面：采用与步骤1a中确定的计算子面大小相同的标定线框，从待测量隧洞衬砌混凝土管片表面的左下角开始，按照从左到右，从下到上的顺序依次进行计算子面的划线标定，并编号，得到m个计算子面；同一行或同一列中计算子面的间隔保持接近或相同。
14.步骤1c、标定特征子面：采用标记笔，从待测量隧洞衬砌混凝土管片表面的左下角开始，按照从左到右，从下到上的顺序依次进行特征子面的轮廓划线标定，并编号，得到n个表征子面；其中，5≤n＜m；同一行或同一列中的特征子面的间隔保持接近或相同；其中，特征子面是指待测量管片表面粗糙起伏度存在肉眼可见差异的区域。
15.步骤2、采集图像：采用拍照设备对完成测量子面标定的待测量隧洞衬砌混凝土管片表面进行拍照，得到待测量管片表面图像；拍照时，待测量隧洞衬砌混凝土管片表面的照度控制在250～400lux；同一隧洞衬砌混凝土管片表面采用同一拍照设备；拍照设备的镜头与隧洞衬砌混凝土管片表面之间的间距控制在30cm～50cm。
16.步骤3、裁剪及预处理图像：将步骤2采集得到的待测量隧洞衬砌混凝土管片表面图像，按照标定的计算子面和特征子面进行逐个裁剪、分类保存，并保留步骤1中的编号；然后，对裁剪后的每幅图像均进行预处理，使其转换为满足设定质量要求的灰度图像；图像预处理完成后，得到n个表征子面图像和m个计算子面图像。
17.步骤4、测量实际粗糙度k：采用粗糙度测量设备，对步骤1标定的n个表征子面的粗糙度分别进行测量，进而得到n个粗糙度测量值k；其中，第i个特征子面的粗糙度测量值为ki，1≤i≤n。
18.步骤5、计算图像特征参数f
p
和fz：针对图像预处理得到的每个表征子面图像和每个计算子面图像，分别计算图像特征参数f
p
和fz；其中，f
p
和fz分别为像素空间灰度值的差均值和均方根。
19.步骤6、选取粗糙度计算通式，具体选取方法如下：
20.步骤6a、当步骤4中n个粗糙度测量值k均不超过2.5mm，且n个粗糙度测量值k相接近或相同时，则粗糙度计算值jrc为关于f
p
的一次线性函数；采用步骤5计算得到的表征子面图像的f
p
和步骤4中对应表征子面的粗糙度测量值k进行拟合，得到一次线性函数的拟合参数。
21.步骤6b、当步骤4中n个粗糙度测量值k为除步骤6a外的其他情形时，粗糙度计算值jrc为关于f
p
和fz的二元多次线性函数；采用步骤5计算得到的表征子面图像的f
p
和fz，以及步骤4中对应表征子面的粗糙度测量值k进行拟合，得到二元多次线性函数的拟合参数。
22.步骤7、计算测量子面粗糙度计算值：根据步骤6选取的粗糙度计算通式，并利用步骤5计算得到的f
p
和fz，对每个计算子面和每个特征子面分别计算粗糙度计算值jrc；其中，
第i个特征子面的粗糙度计算值为jrci；第j个计算子面的粗糙度计算值为jrcj，1≤j≤m。
23.步骤8、确定误差修正系数：将表征子面的粗糙度计算值jrc和粗糙度测量值k进行误差分析，得到误差修正系数ξ，具体计算公式为：
[0024][0025]
步骤9、计算衬砌混凝土管片表面粗糙度具体计算公式为：
[0026][0027]
步骤1a中，计算子面的大小和数量具体根据待测量管片的实际表面面积s进行确定，具体为：
[0028]
当s≤10m2时，选取的所有计算子面的面积之和不得小于s/10。
[0029]
当s＞10m2时，选取的所有计算子面的面积之和不得小于s/100，且同一采集环境下，同一隧洞衬砌混凝土管片表面上的计算子面的数量不得小于10个。
[0030]
步骤1c中，特征子面的数量不得小于计算子面的1/10，且同一采集环境下，同一隧洞衬砌混凝土管片表面上的特征子面的数量为5～20个。
[0031]
步骤2中，拍照设备为具有拍照功能的手机、具有拍照功能的ipad或数码相机。
[0032]
步骤2中，拍照时，当采用人工光源对待测量隧洞衬砌混凝土管片表面进行补光时，人工光源与待测量隧洞衬砌混凝土管片表面之间的夹角为90
±
20
°
。
[0033]
当待测量隧洞衬砌混凝土管片表面的粗糙度超过2.5mm时，人工光源与待测量隧洞衬砌混凝土管片表面之间的夹角为90
±
10
°
。
[0034]
假设步骤3中每个表征子面图像和每个计算子面图像的左下角点为坐标原点，过坐标原点的两条边长分别为x轴和y轴；则步骤5中，f
p
和fz的计算公式分别为：
[0035]fp
＝v
p
/a
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0036][0037][0038]
式中，v
p
为区域d的像素空间体积；d为选定的计算子面或特征子面；a为计算子面或特征子面的面积；f
max
为区域d中点灰度最大值；f(x,y)为计算子面或特征子面上(x，y)坐标点的灰度值；n
x
、ny分别为x轴和y轴上测点的数量；δx、δy分别表示在x轴和y轴上测点的间距；k、l分别为测点在x轴和y轴方向上的次序号，f
k 1,l 1
为第k 1、l 1序号点对应的灰度值；f
k,l 1
为第k、l 1序号点对应的灰度值；f
k 1,l
为第k 1、l序号点对应的灰度值；f
k,l
为第k、l序号点对应的灰度值。
[0039]
步骤6a中，粗糙度计算值jrc为关于f
p
的一次线性函数的表达式为：
[0040]
jrc＝af
p
b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0041]
式中，a和b为一次线性函数的拟合参数。
[0042]
步骤6b中，粗糙度计算值jrc为关于f
p
和fz的二元多次线性函数的表达式为：
[0043][0044]
其中：
[0045]
f(f
p
)＝af
p
b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0046][0047]
式中：a、b、a
′
、b
′
和c均为二元多次线性函数的拟合参数；f(f
p
)为构造粗糙度；g(fz)为起伏粗糙度。
[0048]
步骤6中，一次线性函数的拟合参数和二元多次线性函数的拟合参数在拟合时，相关系数r2需大于0.9、平均相对误差mre小于5％。
[0049]
本发明具有如下有益效果：
[0050]
1、本发明基于图像测距的原理，不会对被测表面造成损伤或者破坏，克服了接触式测量方法的弱点。另外，本发明所采用的测量设备易于获取，操作过程简单，自动化程度高、综合性价比高。
[0051]
2、本发明直接采用像素空间的灰度特征值来描述表面粗糙度，通过分析图像像素空间的灰度值即可测量计算混凝土管片表面粗糙度。计算的粗糙度数值不仅能直观反映混凝土管片表面的粗糙起伏程度，还能通过点云数据将其还原成直观的三维图像，实现三维粗糙度的可视化。另外，本发明通过严格控制图像采集、传输、处理及信息采集过程中的操作步骤和因素水平，既保证数字图像技术在隧洞环境下能够正常使用，获取优质的图像，而且最大程度消减测量过程中的随机误差。
[0052]
3、本发明要求所采集的测量子面随机覆盖于整体表面，采集操作按照一定顺序进行，既保证采集子面能覆盖混凝土管片表面所有的粗糙特征，而且方便之后的复查和补拍，避免了局部空白的出现。
[0053]
4、本发明在粗糙度表征上采取了计算和表征分开并行的模式，增加实地测量局部粗糙度、拟合表征函数的操作，在数字图像技术上结合了其他精度更高，表征更直观的测量方法。通过该步骤，可显著提高数字图像测量技术和被测对象的契合度，从而使测量结果更具有针对性，数值更准确。
[0054]
综上所述，本发明所提出隧洞衬砌混凝土管片表面粗糙度测量计算方法速度快、操作简单、性价比高、针对性强、结果准确度高，具有良好的应用前景和推广价值。
附图说明
[0055]
图1显示了本发明一种隧洞衬砌混凝土管片表面粗糙度的测量方法的流程图。
[0056]
图2显示了本发明中不同拍照设备拍摄图像的灰度值与像素点数量之间的曲线图。
[0057]
图3显示了本发明中不同照度下拍摄图像的灰度值与像素点数量之间的曲线图。
[0058]
图4显示了三种粗糙度平面在不同入射角度的人工光源下拍摄的图像对比图。图4(a)显示了三种粗糙度平面在0
°
时拍摄的图像对比图；图4(b)显示了三种粗糙度平面在90
°
时拍摄的图像对比图；图4(c)显示了三种粗糙度平面在180
°
时拍摄的图像对比图。
[0059]
图5显示了拍摄距离与不同类型灰度值之间的曲线关系图。
[0060]
图6显示了采样间距与计算时间的曲线关系图。
具体实施方式
[0061]
下面结合附图和具体较佳实施方式对本发明作进一步详细的说明。
[0062]
本发明的描述中，需要理解的是，术语“左侧”、“右侧”、“上部”、“下部”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，“第一”、“第二”等并不表示零部件的重要程度，因此不能理解为对本发明的限制。本实施例中采用的具体尺寸只是为了举例说明技术方案，并不限制本发明的保护范围。
[0063]
如图1所示，一种隧洞衬砌混凝土管片表面粗糙度的测量方法，包括如下步骤。
[0064]
步骤1、标定测量子面
[0065]
测量子面包括计算子面和特征子面，具体标定方法包括如下步骤。
[0066]
步骤1a、确定计算子面的大小和数量
[0067]
不同隧洞工程混凝土管片粗糙程度以及表面积存在较大的差异，因此对于具体标定的测量子面的数量没有统一定论，需根据表面粗糙程度是否均匀及所求精度的要求进行合理选择。
[0068]
计算子面为在待测量隧洞衬砌混凝土管片表面均匀分布的子面，也即能均匀覆盖到待测衬砌整体表面，从而降低数据的离散性。
[0069]
本发明中，计算子面的大小和数量具体根据待测量管片的实际表面面积s进行确定，具体为：
[0070]
当s≤10m2时，为保证其基本精度，选取的所有计算子面的面积之和不得小于s/10。
[0071]
当s＞10m2时，可根据实际条件、表面面积以及起伏均匀程度酌情降低子面数量，以节约时间成本，提高测量效率，但选取的所有计算子面的面积之和不得小于s/100，且同一采集环境下，同一隧洞衬砌混凝土管片表面上的计算子面的数量不得小于10个。
[0072]
步骤1b、标定计算子面
[0073]
为方便后期构建的像素面得以和实际的粗糙面相对应，统一测量子面的尺寸和分析精度，因此需要准备标定物提前标定表面。本发明中，采用与步骤1a中确定的计算子面大小相同的标定线框，从待测量隧洞衬砌混凝土管片表面的左下角开始，按照从左到右，从下到上的顺序依次进行计算子面的划线标定，并编号，得到m个计算子面；同一行或同一列中计算子面的间隔保持接近或相同。
[0074]
由于混凝土管片表面常存在弧度，二维图像无法精确记录其表面特征，而且人造光源光线较为集中，照度随着中心会向四周逐渐削减，因此需对标定尺寸做出限定，一般对平整或接近平整的混凝土管片表面，光线较为发散的人造光源优选采取100mm
×
100mm尺寸的细木纯白色制框架，而对于衬砌面曲率在0.003m-1
以上或存在肉眼可见曲度的混凝土管片表面，光线较为集中的则选择50
×
50mm尺寸的细木纯白色制框架。
[0075]
在本实施例中，待测的混凝土管片表面整体面积约为30m2，具体尺寸为3m
×
10m,管片表面为露出骨料的二级配混凝土，表面粗糙程度分布均匀，粗糙度较大。因此选择
100mm
×
100mm的标定线框。确定计算子面为60个，采集顺序为从左到右均匀采集10个子面，从下到上采集6层，依次编号1～60。
[0076]
步骤1c、标定特征子面
[0077]
此外为方便后续图像数字化后提取的特征参数值可以较好的表征实际粗糙度，需要挑选若干表征子面，对这些子面进行实际粗糙度的测量，表征子面选取只需覆盖混凝土管片表面所有粗糙起伏度差异明显的区域(该差异通过目测判别，如存在明显的凹陷区，两块区域的表面粗糙起伏程度存在肉眼可见的差异等)，再进行均匀的选取，特征子面的数量优选不得小于计算子面的1/10，且同一采集环境下，同一隧洞衬砌混凝土管片表面上的特征子面的数量为5～20个。
[0078]
采用标记笔，从待测量隧洞衬砌混凝土管片表面的左下角开始，按照从左到右，从下到上的顺序依次进行特征子面的轮廓划线标定，并编号，得到n个表征子面；其中，5≤n＜m，同时n≤20。同一行或同一列中特征子面的间隔保持接近或相同；其中，特征子面是指待测量管片表面粗糙起伏度存在肉眼可见差异的区域。
[0079]
在本实施例中，表征子面6个，用白色粉笔在壁面勾勒，并编号1-6。
[0080]
步骤2、采集图像
[0081]
采用拍照设备对完成测量子面标定的待测量隧洞衬砌混凝土管片表面进行拍照，得到待测量管片表面图像。
[0082]
上述拍照设备优选为具有拍照功能的手机、具有拍照功能的ipad或数码相机等。手机的型号不限，如华为、小米或苹果等。
[0083]
图2中，(a)图分别显示了camera/500(500是指相机感光度值)、camera/2000、camera/4000和phone四种拍照设备所拍摄图像的灰度值与像素点数量之间的曲线图；(b)图分别显示了camera/500、camera/2000、camera/4000和phone四种拍照设备所拍摄图像的灰度值与像素点数量之间的曲线图。
[0084]
上述图2中的(a)和(b)来自于两个不同的混凝土表面，目的在于验证虽然不同的拍摄设备和内在参数采集到的表面图像灰度直方图的峰值、像素点数量等存在差异，但是其直方图的整体的趋势是一样的，说明对于表面粗糙度的特征描述是相同的。
[0085]
目前已经将直方图换成了图像灰度值分布曲线图。不同的拍照设备，其所成的灰度直方图分布特点均极为接近，趋势基本相似，而其中手机和相机相比，相机的iso值的差异明显，但依旧能很好的表现出粗糙面的特征，因此，每一套设拍照备均可以采集合格的图像数据进行粗糙度的计算；此外考虑到在工程现场测量的随机性即部分条件的约束，本测量方法的图像采集设备可灵活选取，只需要控制同一项测量工作采取同一个设备即可。本实施例中，优选采用型号为eos7d的佳能数码相机。
[0086]
拍照时，光环境和拍摄距离均对最终粗糙度值的测量具有较大影响。其中，光环境包括光的照度和光的入射角度。下面分别对照度、入射角度和拍摄距离的优先值进行介绍如下。
[0087]
1、照度
[0088]
图3显示了本发明中不同照度下拍摄图像的灰度值与像素点数量之间的曲线图。其中，共设置六种类型的光照环境，分别为表示太阳光充裕的正午或午后，光照度在40000lux以上；使用黑暗中相机闪光灯的情形，光照度在40000lux以上；阳光充足的清晨或
下午，光照强度约为9680lux；光线稀薄的阴雨天，光照强度在3240lux；黑暗情况下人造光照，光照度在252lux；光线充足的室内环境，光照度在397lux。
[0089]
当照度在250～400lux之间的人造光和室内光两种方案下，图像所称的各灰度值像素点数量分布和实际粗糙面各坐标点的高度分布规律极为相似，而在3000～10000lux的普通光和较弱光两种方案下，虽大致趋势还比较相似，但已经出现了一定程度的差异，而照度在40000lux以上的剩余两个方案，其直方图分布明显不符合实际。
[0090]
因而，本发明的照度选择为250～400lux，杜绝闪光等强光源下所得图像。该照度人工光源可代替自然光源实现照度要求，此外该照度也符合隧洞环境。本实施例中，照度优选选择为270lux，误差波动约为20lux。
[0091]
2、光照角度
[0092]
由于粗糙表面必然存在着凹凸起伏的微小山丘状，而考虑隧洞环境下，很难做到全方位的均匀光照，从而使粗糙表面会存在大量的局部阴影，干扰粗糙面的图像数据和粗糙面计算，再结合人工光源发射的光源角度可控的条件，可证明研究光线入射角度具有重要意义。本实施例中，人工光源优选为iphone8五级的手机电筒。
[0093]
图4显示了三种粗糙度平面在不同入射角度的人工光源下拍摄的图像对比图。图4(a)显示了三种粗糙度平面在0
°
时拍摄的图像对比图；图4(b)显示了三种粗糙度平面在90
°
时拍摄的图像对比图；图4(c)显示了三种粗糙度平面在180
°
时拍摄的图像对比图。
[0094]
在图4中，第一列为s-1级别的粗糙度平面，第二列为s-2级别的粗糙度平面，第三列为s-3级别的粗糙度平面。粗糙度值排序为：s-1＜s-2＜s-3。从图4可知，垂直90
°
入射时的粗糙度表面最为清晰。
[0095]
因此必须保证入射光线需要垂直照射到粗糙面，但在实际操作中，受条件限制难免会发生偏差，因此对于粗糙度较小的磨光砂面可控制在偏离垂直20度左右；也即人工光源与待测量隧洞衬砌混凝土管片表面之间的夹角为90
±
20
°
。
[0096]
然而，对于粗糙度较大，表面凹凸起伏较为明显的咬合面(如待测量隧洞衬砌混凝土管片表面的粗糙度超过2.5mm)，需要将角度控制在10度以下；人工光源与待测量隧洞衬砌混凝土管片表面之间的夹角为90
±
10
°
。这样所获得的图像才能和垂直照射所获得的图像特征及数据较为接近，从而满足精度要求。
[0097]
本实施例中，光源优选产生光线与测量面垂直。
[0098]
3、拍摄距离
[0099]
拍摄距离在摄影中是一个重点，镜头和实例距离的改变或焦距的变化都会直接影响成像结果。在实际拍摄中，该距离常无法发进行精确控制，总会发生几厘米的波动。
[0100]
图5显示了拍摄距离与各图像特征参数和计算时间之间的曲线关系图。其中，图像参数包括灰度平均值f
mean
、灰度均方根值fz、灰度值标准偏差s
td
、灰度差均值f
p
和参数计算时间time。
[0101]
通过实验发现，当镜头距离粗糙面40cm～50cm时，结果参数和图像的直方图分布变幅较小。因此设定拍摄距离在30cm～50cm之内，本实施例中，拍摄距离优选为40cm。该范围的拍摄距离下，常见的相机镜头均能满足要求，同时可以保证采集的子面图像具有较高的像素水平。光源应该控制在相机镜头前后5cm左右的距离，本实施例中，光源优选位于相机镜头上端约3cm处。这样的间距控制，可以保证相机拍摄时不会受到光源设备的影响，同
时光源发出的光线又可以近乎垂直照射于混凝土管片表面。
[0102]
这样所得图像的质量不会因为拍摄过程中拍摄距离的变化发生较大程度的偏差，且该距离下市面上流通的相机镜头都能直接实现清晰的摄像。而距离若减小，有些镜头便无法成像，拉大则会导致图像模糊，质量直线下降。
[0103]
步骤3、裁剪及预处理图像
[0104]
运行matlab程序，将步骤2采集得到的待测量隧洞衬砌混凝土管片表面图像，按照标定的计算子面和特征子面进行逐个裁剪、并按子面的类型分开保存、操作，同时保留步骤1中的编号。
[0105]
然后，对裁剪后的每幅图像均进行预处理，使其转换为满足设定质量要求的灰度图像。
[0106]
研究发现常规拍摄所得的子面图像一般为rgb三通道的图像，虽颜色鲜艳，但若想用图像数据直接描述混凝土管片表面的粗糙度，单通道的灰度值较之三通道的rgb数值更加简洁明了；此外在采集和传输过程中由于一些不可控因素，图像会出现不同程度的噪声、局部发虚等问题，若不做处理便直接提取图像数据，结果会产生较大的误差，因此在数据提取之前，需对所采集的图像进行预处理。
[0107]
本发明的图像预处理主要包括图像灰度化、图像降噪以及锐化处理。图像灰度化采取均值化理论，处理过程通过编写matlab程序来实现。
[0108]
平滑降噪处理优选高斯滤波。高斯降噪的基本原理为扫描图像的每一个像素点，以其值及周边领域内其他像素点的值加权平均后代替原值。通过高斯滤波主要可以消除高斯噪声，该噪声由电路元件、图像传感器的自身缺陷导致，为图像噪声的主要组成部分，消除高斯噪声可以大大降低噪声对图像信息的干扰。
[0109]
锐化处理优选sobel算子。混凝土管片表面起伏的凸起和下凹边缘对于描述其粗糙程度极为重要，成像时，其往往成为像素点数值发生跃迁的地方。因此，为避免这种边缘特征被模糊化，需要采取锐化处理。sobel算子根据像素点上下、左右邻点灰度值的加权差在边缘处达到极值这一现象检测边缘，而且对噪声亦能起到平滑作用。
[0110]
本发明采用gauss-sobel的组合算法显著降低了图像的输入和传输误差，两种算子相辅相成，使图像数据更加真实的反映表面粗糙度。由于混凝土管片表面粗糙度、测量的环境为较为嘈杂的施工环境，其图像噪声较之其他工作环境更为严重，而且工程应用中精度控制一般为毫米级，因此gauss-sobel的组合算法可为优选方案。
[0111]
在进入图像信息的提取前，为避免操作不当产生劣质图像对最终测量结果造成误差干扰，故须进行2道图像筛选：一道于预处理之前，对存在明显缺陷的图像进行删除、复核、复拍等操作；第二道在预处理之后，对预处理结束后的图像进行检查，剔除预处理效果不良的图像。
[0112]
基于以上分析，隧洞衬砌混凝土管片表面图像预处理及筛选的主要操作步骤优选方案为：
[0113]
对裁剪产生的图像进行人工复查，对不合格的图像进行删除，并记录该标号。
[0114]
运行程序，对图像进行图像预处理。
[0115]
对产生在一个新文件内的所有处理后灰度图像进行审查，对不合格的进行删除，并记录编号；综合编号上的所有问题图像，结合采集时的分布情况综合分析，若问题图像数
量较大，或者消除这些图像后，对混凝土管片表面的采集出现局部区域空白，则实行逐级往前复查并重新操作；若未出现以上两个现象则可直接剔除，不再重复操作；确认所有预处理好的图像安置于一个文件夹，保留各自初始编号。
[0116]
需注意的是，所有图像的删除操作仅针对当前检查的文件夹，之前的图片信息保留，此外，若复查后重复操作裁剪、预处理等步骤，则直接覆盖图像，不另起编号。
[0117]
对裁剪后的图像进行预处理、数字化、和像素数据提取，期间经过二道图像筛选。编号18、31的计算子面图像质量及数据不符合要求，表征子面均满足要求，删除18、31的图像并不会降低测量结果的精度，因此不做补拍处理。根据需要，也可进行补拍处理。
[0118]
另外，对比图像的3d扫描所获得的点云数据，发现对于50
×
50mm的粗糙面扫面后所形成的数字面其点云数据量在17000条数据左右，而图像法较之扫描法抗干扰的能力较差，因此适当提高点云数据的数量，100
×
100mm或50
×
50mm的粗糙面图像用来计算分析的数据量不得小于80000/20000条，否则视为劣质数据剔除。本实施例中，数据文本的数据量均在8000条之上，满足数据分析的要求。
[0119]
然后，采用sufer 15对图像数据文本继续插值分析，采样间距选择0.1mm，插值方法选择最小曲率法，单次分析时间约为0.5s，所有图像数据的分析操作时间约为300s，生成数据报告保存为txt格式。
[0120]
上述sufer在对数据进行插值计算时候，采样间距会影响到计算的时间和计算的结果，在图6中，其显示了采样间距与计算时间的曲线关系图。由图6可知，采样间距越小，精度越高，但时间成本也越昂贵。研究发现对于100mm
×
100mm粗糙面图像采样间距选择0.1mm性价比最高，采样间距进一步提高其精度的提升较为轻微但时间确呈现指数型增长，反之，则相反。
[0121]
sufer15的插值方法有克里金插值法、最小曲率法、改进谢别德法、自然邻点插值法、和径向函数插值法等十余种方法，我们对其中经典的九种分析方向，其计算结构均较为接近，偏差极小，但时间差异确较大，在几秒到几十秒不等，因此我们放弃选择传统的克里金插值法(43s)而选择最小曲率法(0.5s)。
[0122]
图像预处理完成后，得到n个表征子面图像和m个计算子面图像。
[0123]
假设步每个表征子面图像和每个计算子面图像的左下角点为坐标原点，过坐标原点的两条边长分别为x轴和y轴。
[0124]
步骤4、测量实际粗糙度k：采用粗糙度测量设备，对步骤1标定的n个表征子面的粗糙度分别进行测量，进而得到n个粗糙度测量值k；其中，第i个特征子面的粗糙度测量值为ki，1≤i≤n。
[0125]
上述粗糙度测量设备，采用目前行业公认的高精度粗糙度测量装置和技术，包括但不限于千分尺测量、3d激光扫描等。本实施例中，优选选择型号为0-1mm/0.001mm的千分表。
[0126]
步骤5、计算图像特征参数f
p
和fz[0127]
针对图像预处理得到的每个表征子面图像和每个计算子面图像，分别计算图像特征参数c和fz；其中，f
p
和fz分别为像素空间灰度值的差均值和均方根。
[0128]
通过sufer可算出包括最大值，最小值，平均值，均方根，差均值，标准偏差等众多统计学参数数值，这些数值可以最直观的反映图像像素空间中各像素点的灰度值大小，通
过相关的函数计算可以来表征粗糙面粗糙程度。
[0129]
采用matlab软件获取图像空间坐标(x，y)对应的像素值f(x，y)，得到图像空间曲面像素模型，即
[0130]
f＝f(x,y)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0131]
式中：x、y分别为粗糙面上测量点的坐标值，mm；f或f(x,y)为计算子面或特征子面上(x，y)坐标点的灰度值。
[0132]
利用式(1)求出粗糙面上给定区域d内各点的像素空间曲面像素模型，进而可计算灰度极大值f
max
所在平面间所围体积v
p
，即
[0133][0134]
式中：v
p
为区域d的像素空间体积；d为选定的计算子面或特征子面；a为计算子面或特征子面的面积；f
max
为区域d中点灰度最大值。
[0135]
图像灰度差平均值为像素空间体积与给定区域d的面积a之比，即
[0136]fp
＝v
p
/a
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0137]
灰度均值f
mean
、灰度起伏度均方根fz及灰度差平均值f
p
的关系为：
[0138][0139]fp
＝f
max-f
mean
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0140][0141]
式中，n
x
、ny分别为x轴和y轴上测点的数量；δx、δy分别表示在x轴和y轴上测点的间距；k、l分别为测点在x轴和y轴方向上的次序号，f
k 1,l 1
为第k 1、l 1序号点对应的灰度值；f
k,l 1
为第k、l 1序号点对应的灰度值；f
k 1,l
为第k 1、l序号点对应的灰度值；f
k,l
为第k、l序号点对应的灰度值。
[0142]
上述每个测点其实就是一个坐标点，用于灰度值读取。测点选取和提取，均通过sufer软件自动完成，仅需要输入间距和插值方法。
[0143]
下面以100mm
×
100mm图像为例，对测点的选取进行详细说明如下。
[0144]
输入间距为1mm，则将图像四条边(也即四个坐标轴)均分为100等分，从而出现平行于x轴与y轴的九十九条线，这些线和坐标轴的交点便是测点，其各自相交便是处于内部的序号点，总共就产生100
×
100共10000个点。在测点选取时，间距越小，数量越多，因而计算的结果也更加精确，但计算的时间成本也越高。本发明中，采取sufer软件对图像数据进行处理，采用0.1mm的采样间距和最小曲率的插值方法，计算效率最高，精度也高。
[0145]
本发明中，在表征通式选取之前的各项工作和设定，目的有如下三点：
[0146]
1、建立起图像参数和粗糙面几何参数之间的良好的映射关系。
[0147]
2、降低操作过程中随机误差的产生，从而为最后统一的误差修复提供一种可能。
[0148]
3、在满足前两点目标的前提下进一步的降低了测算成本，为工程中大范围的测算应用提供一种选择。
[0149]
另外，本发明的图像预处理后，为灰度图，通过绘制灰度直方图，可以看到每一个灰度值(0～255)的像素点数量分布情况，对应的是实际粗糙面最大下凹度均等分256份，得
到256个下凹度数值，而每一个下凹度的坐标点的数量分布情况。因此通过分析图像的灰度值直方图就可以直观、定性的判定在这个图像采集的环境、操作下所得到图像是否能反映实际粗糙面的粗糙程度和通过数值计算得到的图像结果参数能否用来表征实际粗糙度，这个表征误差是多少。
[0150]
步骤6、选取粗糙度计算通式
[0151]
由之前的研究发现，目前关于研究表面粗糙度标准的公式有很多，但是其适用性都不太广，总体上分为三大类。构造深度、曲线分形维数法和曲面分形维数法。三种方法都是直接对粗糙面实体或者扫描面进行分析计算，具有较高的精度和实用性，但延用到数字图像中则需要进行对比修正。
[0152]
这个修正后的式子则直接和测量环境和粗糙面的特征有关，不同的相机或手机其所得到的图像特征并不是完全一致的，若直接沿用其他工程或研究论文上的公式难免会水土不服，测量误差增大。对于不同批次、粗糙程度或图像采集环境的衬砌表面粗糙度采用同一种公式去表征，必然会增大系统误差，因此最合适的方法就是采取既有的几种概念，从常见的几种通式选取一种进行拟合。
[0153]
本发明中，粗糙度计算通式的具体选取方法如下：
[0154]
步骤6a、当步骤4中n个粗糙度测量值k均不超过2.5mm，且n个粗糙度测量值k相接近或相同时，则粗糙度计算值jrc为关于f
p
的一次线性函数；采用步骤5计算得到的表征子面图像的f
p
和步骤4中对应表征子面的粗糙度测量值k进行拟合，得到一次线性函数的拟合参数；其中，n个粗糙度测量值k相接近是指任意两个粗糙度测量值k之间的差值不超过两者中较大值的1/2。
[0155]
上述粗糙度计算值jrc为关于f
p
的一次线性函数的表达式为：
[0156]
jrc＝af
p
b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0157]
式中，a和b为一次线性函数的拟合参数，在拟合时，相关系数r2需大于0.9、平均相对误差mre小于5％。
[0158]
步骤6b、当步骤4中n个粗糙度测量值k为除步骤6a外的其他情形时，粗糙度计算值jrc为关于f
p
和fz的二元多次线性函数；采用步骤5计算得到的表征子面图像的f
p
和fz，以及步骤4中对应表征子面的粗糙度测量值k进行拟合，得到二元多次线性函数的拟合参数。
[0159]
上述粗糙度计算值jrc为关于f
p
和fz的二元多次线性函数的表达式为：
[0160][0161]
其中：
[0162]
f(f
p
)＝af
p
b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0163][0164]
式中：f(f
p
)为构造粗糙度；g(fz)为起伏粗糙度；a、b、a
′
、b
′
和c均为二元多次线性函数的拟合参数；在拟合时，相关系数r2需大于0.9、平均相对误差mre小于5％。
[0165]
当粗糙度计算值采用二元多次线性函数，也即公式(10)进行计算时，其中，起伏粗糙度g(fz)计算时，需要对公式(9)中的三种通式进行依次计算，然后选取通式的拟合结果
最好，即mre最小的通式，作为公式(10)所需要的起伏粗糙度值。
[0166]
本实施例中，优选步骤6a中的一次线性函数，图像表征参数选为像素灰度差均值，提取表征子面图像的像素灰度差均值和实际粗糙度进行线性函数的拟合，r2为0.98，满足拟合精度的要求。
[0167]
步骤7、计算测量子面粗糙度计算值：根据步骤6选取的粗糙度计算通式，并利用步骤5计算得到的f
p
和fz，对每个计算子面和每个特征子面分别计算粗糙度计算值jrc；其中，第i个特征子面的粗糙度计算值为jrci；第j个计算子面的粗糙度计算值为jrcj，1≤j≤m。
[0168]
步骤8、确定误差修正系数：将表征子面的粗糙度计算值jrc和粗糙度测量值k进行误差分析，得到误差修正系数ξ，具体计算公式为：
[0169][0170]
本实施例中，使用上述公式(11)，确定误差修正系数为0.9。
[0171]
步骤9、计算衬砌混凝土管片表面粗糙度具体计算公式为：
[0172][0173]
以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。