杨 阳，张 亮，张洪军，王 池，孔海燕，朱泽军

(1.中国计量大学 计量测试工程学院，浙江杭州310018；2.中国计量科学研究院，北京100029； 3.郑州计量先进技术研究院，河南郑州450001；4.河南省环境监测技术重点实验室， 河南郑州450004；5.河南省生态环境监测中心，河南郑州450004)

1 引 言

我国采取了一系列行动来治理大气污染，建立并完善区域大气污染联防联控机制，在《中华人民共和国大气污染防治法》中指出：固定源烟道排放是大气污染物的最主要来源[1～4]。

为了实现对烟道大气污染物排放量的监管与控制，对于大中型固定排放源通常在烟道上安装连续排放监测系统(continuous emission monitoring system,CEMS)[5,6]，其中包括使用烟气流量计监测烟气流量。但存在无法直接将现场烟气流量计溯源到现有实验室流量计量标准装置的问题。针对上述情况，本课题组研制了基于速度面积法的现场烟气流量在线校准装置，实现对安装在烟道上的流量计进行在线校准[7,8]。

速度面积法是指在一个烟道截面上测量多个速度点，通过速度积分计算截面流量的方法。目前，国际上基于速度面积法的烟道流量在线比对方法有：EPA(environmental protection agency )方法,该方法提出将烟道截面划分若干相等区间，将测量点置于这些相等面积内，通过计算测量点流速平均值获得截面平均流速; 还有在ISO 7194和ISO 3966标准中介绍的一种使用皮托管对封闭管段流体流动的基于速度面积法的测量方法；在GB 16157中规定的固定污染源烟道流量现场比对方法，该方法类似EPA方法。

通常管道横截面积数据获得方式包括：通过设计图纸获得；通过激光测距仪实际测量烟道内部特征尺寸，使用全站仪或激光扫描仪测量烟道外部尺寸后根据壁厚计算内部尺寸等获得。前者由于实际烟道尺寸和设计图纸的差异以及长时间管道内壁结垢等原因，实际烟道横截面积与图纸尺寸会有出入，且这一差异难以估计。后者测量时激光测距仪的激光通过法兰口打入，通过测量由管道内壁反射回的激光光束得到管道直径，由于手持激光测距仪难以保证打入方向为直径方向，角度偏差会对测试结果产生较大影响，且对于矩形烟道往往只能获得截面的一个边长，无法计算截面面积；如使用全站仪和激光扫描仪进行测量，测试时企业很难配合停机、清理烟道并将仪器放入管道内进行测量，因此往往测得管道外尺寸，再通过测厚仪测量管道厚度，从而计算烟道内部截面面积，此方法流程复杂且烟道壁厚测量可能不均匀，因而测量准确度较低。

利用烟道截面校准装置能较好解决上述问题，装置由法兰口进入烟道后，激光探头会按设定角度旋转，其中包括平行于横截面方向的360°旋转和垂直于横截面的上下摆动，从而得到激光探头距离整个烟道内壁的距离和角度，可转换为烟道内壁上各点相对旋转轴的球坐标三维点云，再通过选择不同的数学模型拟合得到圆柱方程，从而求得截面面积。圆柱面拟合方法主要包括遗传算法、随机抽样一致算法、最小二乘法、特征值算法、距离函数参数化算法等[9]。本研究主要比对了随机抽样一致性圆柱拟合算法[10,11]和最小二乘法圆柱拟合算法[12,13]，获得最优模型，评估了烟道截面校准装置的测量偏差和不确定度[14,15]，结果令人满意。

2 校准装置介绍和测量原理

2.1 装置介绍

图1为烟道截面面积校准装置的示意图。装置主要由激光探头，旋转机构，弯折机构和控制箱组成。装置可以插入烟道法兰孔并固定安装，测量前使用弯折机构将装置前半部分旋转90°至俯视图所示位置。激光探头为高精度的激光测距传感器，旋转机构能够改变激光探头的俯仰角(-60°～60°)和扭转角(360°任意旋转)，实现对烟道的完整测量。控制箱则保证了激光探头和旋转机构的自动化测量。

图1 烟道截面面积校准装置示意图Fig.1 Schematic diagram of stack cross section area calibration facility

2.2 测量原理

图2是烟道截面面积校准装置在DN2000管道的测量照片。量臂由法兰孔插入烟气管道内，使用大力钳将其稳固并锁定。使用弯折机构将激光探头旋转90°至图1所示俯视图位置，使用控制软件控制系统进行逐点扫描测量，测量数据为以旋转轴轴心为原点的球坐标点云。构建的球坐标系中以旋转轴轴心为原点，激光探头摆动方向为俯仰角，旋转轴旋转方向为方向角。测量时，装置会自动调整到起始俯仰角，并在此俯仰角下做方向角的360°逐点扫描，随后再逐渐改变俯仰角完成所有测试，测量中激光探头到烟道壁面各点距离需要加上激光探头与旋转轴轴心的距离，方可得到坐标原点到烟道内壁各点的距离。

图2 烟道截面校准装置测试Fig.2 Stack cross section area calibration facility test

为了方便数据处理再将球坐标转为笛卡尔坐标，转换后得到烟道内壁三维点云的原始数据，图3为测量数据的空间成像。

图3 测量数据的空间成像Fig.3 Spatial imaging of measurement data

3 烟道截面面积拟合模型选择

3.1 随机抽样一致性算法

在模型估计领域，算法有着十分重要的地位，随机抽样一致性算法是通过随机抽样与内点计算而确定最优模型的最优化估计算法[16]。同时为了找到最优模型，适当地割舍掉了一些现有的数据。

3.1.1 随机抽样算法原理和特点

随机抽样一致性算法的原理如下：首先在已经得到的原始数据中随机选择几个点作为内群点；根据内群点的数学关系求解模型方程；再将所有的样本点代入模型方程中算出各点误差；设置误差上下限并记录在此误差限以内的样本点个数记为内群点的个数；重复随机抽样过程并选择拟合效果中内群点个数最多的方程为最终拟合方程。

该算法根据其拟合原理具有以下特点：

(1)拟合效果极大地避免了干扰数据的影响并能根据误差方程反馈出内群点的个数；

(2)拟合效果与随机抽样重复次数直接相关；

(3)拟合效果存在一定的不确定性。

3.1.2 随机抽样在圆柱拟合的应用

随机抽样拟合中需要求得原始数据各点的单位法向量，首先搜索如图3所示任意一点的若干临近点，再将这些点拟合成一个平面，所得到的该平面的法向量即为该点的单位法向量。图4为计算得到的原始数据各点的单位法向量及三维点云重建。

图4 三维点云重建Fig.4 Three-dimensional point cloud reconstruction

然后从包含有单位法向量的原始数据中，随机抽取两个点，根据两点及其单位法向量可以拟合求得圆柱方程；再将原始点三维坐标代入圆柱方程中求出各点与轴线的距离；根据管道半径R的1%误差设置误差限(由于随机抽样算法中忽略了部分数据，因此不能选用总数的3倍标准差)，找出在此误差限内的内群点的个数；多次重复随机抽样流程，找出内群点最多的一组，确定为圆柱拟合最终方程。

3.1.3 随机抽样拟合效果对比

根据随机抽样一致性算法的特点，对比不同随机抽样次数N对拟合效果的影响，结果如图5所示。

图5 抽样次数对拟合效果的影响Fig.5 The influence of sampling times on the fitting effect

3.2 最小二乘法

3.2.1 剔除离群值及比对

烟道截面面积基准装置在实际测量过程中因为激光探头旋转时光斑打到测量臂上或者通过管道法兰口打出管道，出现了一部分异常数据。为不影响后续工作的数学处理，根据推算，将误差超过测量距离平均值0.5倍的数据予以删除，图6为剔除前后成像对比。

图6 剔除粗大误差数据前后效果对比Fig.6 Comparison before and after removing gross error data

3.2.2 最小二乘法的圆柱拟合

最小二乘法圆柱拟合的原理是通过求解圆柱的轴线方程，再将数据点按轴线方向投影得到一个在投影面的圆，最后对圆进行最小二乘曲线拟合从而得到投影圆的半径，回带求出圆柱方程完整信息。具体流程如下：

分别求出X,Y,Z三轴的重心坐标，由于圆柱体的重心一定在圆柱体的轴线上，可以建立特征方程，求解特征方程可以得到圆柱体的轴线方程；将之前求得的笛卡尔坐标系进行坐标转换，使得Z轴和圆柱的轴线方程重合；坐标轴按Z轴投影即可以得到近似一个圆；通过最小二乘法曲面拟合，可以求得圆的圆心及半径；最后将圆心坐标进行坐标反变换求得圆柱轴线上一点，根据圆柱轴线、轴线上一点和半径可以构建圆柱标准方程。

3.2.3 迭代剔除误差并比对

虽然在第3.2.1节已经将初步推算的物理遮挡或测量到非柱面点造成的离群数据部分删除，但是仍然有非壁面点的数据残余，因此提出一种迭代的思想剔除干扰数据。

基于3倍标准差对数据进行筛选剔除并重新拟合，当原始数据点个数为11 041时，如上图7所示拟合效果在第5次迭代时收敛并趋于稳定，当原始数据点个数增加时，拟合效果会在更高的迭代次数下收敛并趋于稳定。

图7 迭代次数对拟合效果的影响Fig.7 The effect of iteration times on the fitting effect

4 校准装置测量不确定度评定

通过对比上述两种算法拟合效果，最终选择最小二乘法圆柱拟合，并对基于此方法的烟道截面面积校准装置测量不确定度进行评定，评定标准参照国标GB 34881(坐标测量机的检测不确定度评估指南)。装置测量不确定评定主要包括以下4部分：烟道截面面积校准装置的A类测量不确定度；激光探头的测距不确定度；单一点测量稳定性的不确定度；温度变化导致管道半径变化的不确定度。

4.1 部分测量不确定度

(1)截面面积校准装置的A类测量不确定度

使用烟道截面校准装置对DN2000管道半径进行3次测量，得到R1=1 005.21 mm,R2=1 005.42 mm,R3=1 005.37 mm。根据测量数据能够计算得到烟道截面面积校准装置A类测量不确定度，见表1。

表1 烟道截面面积校准装置测量数据Tab.1 Measurement data of stack cross section area calibration facility

由表1 中的3次测量数据可得出，相对扩展不确定度Ud=0.013%，k=2。

(2)激光探头的测距不确定度

激光探头安装于装置前端，用来发射激光光束，并对测量数据进行采集和记录，其有效测量范围是1～20 m。装置经中国计量科学研究院长度所校准,激光探头测距的扩展不确定度为U1=0.7 mm(k=2)，标准不确定度为u1=0.35 mm。测量过程中的管道半径取R=1 005.33 mm,则可计算激光探头的相对标准不确定度为uk=0.035%。相对扩展不确定度Uk=0.070%，k=2。

(3)单一点测量稳定性的不确定度

单一点测量稳定性的不确定度，是对仪器的稳定性考核，考察测量同一位置的时候，数据的漂移量，以及数据稳定性随时间的关系，图8是持续间隔20 min对同一位置多次测量的结果。

图8 时间变化对测量结果的影响Fig.8 The effect of time change on the measurement results

(4)温度变化导致管道半径变化引入的不确定度

测试过程中实验管道为碳钢材料，需要考虑其热膨胀系数α的不确定度。温度变化导致管道半径变化的标准不确定度为

uc=R·α·ut

(1)

式中：R=1 005.33 mm为测量过程中的管道半径平均值；α=0.000 017为碳钢热膨胀系数;ut为测量过程中温度引起的不确定度，ut=Δt/k,其中Δt取多次测试中最大温度变化为3℃，扩展因子k=3;因此，得到：uc=0.017 mm。相对扩展不确定度Uc=0.005 1%，k=3。

4.2 装置测量截面面积不确定度合成

4.2.1 装置测量半径不确定度合成

装置测量半径不确定合成主要包括: A类不确定度烟道截面面积校准装置的测量不确定度ud;B类激光探头的测距不确定度uk，单一点测量稳定性的不确定度uw，温度变化导致管道半径变化引入的不确定度uc，其合成不确定度为：

(2)

取k=2,相对扩展不确定度UR=0.084%。

4.2.2 装置测量截面面积不确定度合成

装置测量截面面积的不确定度需要从半径合成到截面圆面积，其合成关系如下：

uS=2uR

(3)

装置测量面积相对标准不确定度uS=0.084%，相对扩展不确定度US=0.17%，k=2。

5 结 论

烟道截面面积校准装置是烟道流量在线校准装置的重要组成部分，保证了基于速度面积法烟道流量在线校准中烟道横截面积的准确度和可靠性。本文对烟道截面面积校准装置面积拟合模型的特点和拟合效果进行深入研究，对装置的测量不确定度进行评定，得出以下结论:

(1)拟合模型中随机抽样一致性算法的拟合效果随着随机抽样点次数增加而变得更好，该方法能避免误差数据的干扰，当抽样点次数为样本10倍以上时，效果趋于稳定，此方法拟合的平均值最大偏差为0.27%。

(2)拟合模型中最小二乘算法的拟合效果应避免粗大误差数据的干扰，采用离群点删除和3倍标准偏差迭代剔除错误数据，拟合效果随迭代次数的增加收敛并趋于稳定，方法拟合结果的平均值最大偏差为0.02%。

(3)选择最小二乘法为烟道截面标准装置的拟合方法，优化后的数据测量重复性为0.013%。同时对烟道截面标准的不确定度进行评定，结合其它因素合成装置截面面积测量不确定度为0.17%。