Locating algorithm of optical fiber spot center based on FPGA

随着科学技术的进步，光纤被广泛应用在各个领域，光纤光斑中心位置检测是光电测量手段中的关键技术，其准确性、精度和抗干扰性直接影响测量结果。目前光斑中心定位算法主要有Hough变换算法、重心法和拟合法。鉴于Hough变换算法需要逐点进行投票、记录，花费时间较长；重心法则要求光斑图像光强分布均匀，且其受光斑形状的影响较大，会产生误差；而拟合法时间复杂度较小、运算精度高^[1]，但容易受噪声影响。

随着集成电路的发展，可编程门阵列(field-programmable gate array, FPGA)器件有着以下特点：优秀的可编程能力、并行执行能力、硬件资源丰富、低功耗。因此FPGA被广泛应用于图像处理算法的领域。现有测量激光光束质量分析的系统都已经相当成熟，但都是通过计算机对激光光斑图像进行处理来计算光斑的参量。例如，参考文献[2]中基于LabVIEW平台开发了一套快速检测激光发散角的测量系统。而参考文献[3]中将FPGA中采集到的激光光斑数据经USB线与上位机通信，并对数据二次处理，得到光斑的能量分布图，并解析出光斑的中心。参考文献[4]中对Hough变换法、重心法以及高斯拟合算法定位光斑中心相比较，得到高斯拟合方法效果最好，误差小于0.1pixel。参考文献[5]中对传统最小二乘法拟合圆以及改进的最小二乘法拟合圆进行比较，改进的圆拟合方法处理得到的中心坐标误差在0.1pixel之内。因此，为保证误差小于0.1pixel并能够提高计算效率，研究基于FPGA平台处理光斑数据的系统是非常必要的。针对这一问题，相对于传统软件进行光斑图像处理的算法，本文中提出一种基于FPGA的光纤光斑中心定位的算法，提高了运算速率，并且增加了图像预处理部分，可高效地过滤图像中的噪声，采用Sobel算子提取光斑轮廓，获得光斑边界像素坐标，再利用代数距离的最小二乘逼近光斑的轮廓，从而得到光斑的中心位置坐标，为求取光斑质量参量奠定基础。

目前，常见的中心定位方法有Hough变换、重心法和最小二乘拟合^[6]。鉴于重心法和Hough变换法具有复杂度高和精度低的缺陷，因此，本文中选用最小二乘法拟合光斑边界。最小二乘法拟合算法对噪声敏感，所以，问题的根本在于提高该算法的抗干扰能力，以提高光斑中心定位的精度。

中值滤波是一种非线性平滑技术，通常用于处理椒盐噪声。其实现过程为：首先对图像中的某个采样窗口中的所有数据进行排序，然后用排序后的中值作为该区域的中心位置的像素值。传统的中值滤波算法运算复杂度高，故本文中采用快速中值滤波算法：在3×3的滑动窗口中，先比较每列值中的最大值、最小值、中值；再取出3个最大值中的最小值A₀，3个中值的中值A₁，3个最小值中的最大值A₂；最后取出A₀, A₁, A₂的中值即为最终结果。与线性滤波相比，中值滤波可以保护图像的边缘信息并且有效滤除噪声。经滤波以后的图像如图 1所示。

Figure 1.  Spot image before and after filtering

图像经中值滤波后，需对图像进行阈值分割，即二值化处理——选取适当的阈值^[7]，用0或255将灰度图像分成黑与白两种颜色，本文中选取的阈值为图像中灰度最大值的一半，而后采用边缘检测算法提取光斑的轮廓。Sobel算子与Prewitt算子相比，对像素位置的影响做了加权，降低了噪声的影响，用其进行目标光斑的边缘检测。换言之，同时测量水平方向和竖直方向的灰度差，将得到的数据叠加。经Sobel检测后，会得到所需的边缘像素坐标，再对数据进行处理。经过Sobel算子边缘检测后的图像如图 2所示。

Figure 2.  Spot image before and after edge detection with Sobel operator

对图像数据进行预处理后，得到平滑的光斑边缘，进而采用曲线拟合对光斑的中心精确定位。圆可以看成是一个特殊的椭圆，本文中研究一般椭圆的光斑中心定位^[8-9]。由图像预处理得到的所有的边缘数据对圆进行最小二乘拟合，利用求解出的方程系数可直接计算出光斑的中心位置。

在2维坐标系中，平面椭圆方程的隐式方程表示为F(a, x)=0，其中$a = (A, B, C, D, E, F) $为椭圆方程的参量，即椭圆隐式方程为：

式中，$ f\left( {x, y} \right)$为椭圆方程，其中x和y分别为椭圆上点的横纵坐标值; $ A, B, C, D, E, F$为椭圆方程参量。

比较常见的椭圆拟合方法有基于代数距离的最小二乘法。由最小二乘法可知，椭圆的拟合可以近似看成是求光斑边缘上的i个数据点$\left( {{x_i}, {y_i}} \right) $到曲线$ f\left( {x, y} \right) = 0$的残差平方和的最小值，即求目标函数$f\left( {B, C, D, E, F} \right) = \sum\limits_{i = 1}^m {} {\left( {{x_i}^2 + B{x_i}{y_i} + C{y_i}^2 + D{x_i} + E{y_i} + F} \right)^2} $的最小值。再由极值原理，欲使$f\left( {B, C, D, E, F} \right) $值为最小，必有$ \frac{{\partial f}}{{\partial B}} = \frac{{\partial f}}{{\partial C}} = \frac{{\partial f}}{{\partial D}} = \frac{{\partial f}}{{\partial E}} = \frac{{\partial f}}{{\partial F}} = 0$。由此可得一个线性方程组，将其转换成求线性矩阵方程Ax=b的解。这里A为5×5维实矩阵，x为待求5维列向量，b为n维实向量。

式中，$ \left( {{x_i}, {y_i}} \right)$为光斑边界上点的坐标，$\left( {{x_i}, {y_i}} \right) \subset E $表示光斑图像中边缘上所有的点的集合。通过求解该线性方程组，便可得到椭圆方程参量B, C, D, E和F，则中心坐标$\left( {{x_{\rm{c}}}, {y_{\rm{c}}}} \right) $为：

在数字信号处理(digital signal processing，DSP)系统设计开发工具中以高斯消元法为基础对线性方程组进行求解，鉴于DSP builder的除法模块只能计算出被除数和除数相除以后的商和余数，故而采用移位寄存器将被除数移位倍增，再进行相除，经验证，其精度远远小于单精度浮点数，计算的结果与理论比较误差大于10%，因此选用Quartus Ⅱ的单精度浮点运算知识产权(intellectual property，IP)核，运用IP核设计运算处理单元。

一般地，线性最小二乘问题的求解主要有迭代法和矩阵分解法两类。其中矩阵分解法可以最大限度实现并行计算，从而更适用于FPGA实现^[10]。矩阵分解的方法主要有QR分解^[11-12]、LU分解^[13]、奇异值分解^[14]以及Cholesky分解^[15]。QR分解的计算过程复杂；而奇异值分解与QR分解的计算难度相当；大规模矩阵一般采用LU分解来进行分块分解，但是LU分解过程中需要大量的乘法和除法运算，会造成一定的延时，并且消耗过多的硬件资源；LDL^T分解适用于对称的正定矩阵，分解过程中避免了开方运算，其计算的复杂度仅为LU分解的一半，故而本文中采取LDL^T分解算法来实现矩阵的分解。

LDL^T分解适用于对称的正定矩阵，对任意一个对称正定矩阵A，将其分解成A=LDL^T的形式，其中L为对角线为1的下三角矩阵，D为对角阵，L^T为L的转置，即：

矩阵D及L中的元素d_j和l_ij根据下式计算:

式中，$ j = 1, 2, \ldots , n;i = j + 1, j + 2, \ldots , n$。a_ij和a_jj分别为矩阵A中的第i行j列和第j行j列所对应的元素，d_j和d_k为对角阵D中对角线上所对应的的元素(k=j-1)，l_ij，l_ik和l_jk分别为矩阵L中第i行j列、第i行k列和第j行k列所对应的元素。

分解后即将求解的线性方程组Ax=b转换为求LDL^Tx=b的解，等同于求解3个三角线性方程组。

式中, z和r为引入的n维待求列向量。根据以上步骤即可求得方程的解向量x。

为了实现在FPGA硬件平台上高效的矩阵分解，必须合理采用流水线技术，设计并行计算^[16]。算法并行化的根本是数据之间的依赖关系，下面将分析LDL^T分解算法的数据之间的关系。

根据(5)式，每一个内层循环用于计算矩阵L的一列元素和矩阵D对应的一个元素，以第j个内层循环中计算的具体数据依赖关系举例说明。将(5)式展开如下所示：

由(7)式可知，矩阵D中元素d_j的计算依赖于矩阵A中的元素a_jj、矩阵L的第j行元素以及矩阵D的前(j-1)个元素。L矩阵的元素l_ij的计算依赖于矩阵A中的元素a_ij、矩阵D中的元素d_j以及矩阵L的第i行元素。而l_ij与l_(i+1)j之间的计算没有依赖关系，且d_j和l_ij的计算都要计算中间结果l_jkd_k。

通过分析数据依赖关系可知，LDL^T分解算法的外层数据循环计算需要依存已经获得的数据，内层数据循环计算则没有依赖关系，直接并行计算。

通过对LDL^T分解算法及相关数据分析，可设计如图 1所示的分解结构框图。分解模块主要分为计算对角阵D元素的单元PED、下三角阵L元素的单元PEL和数据选择器，PED与PEL共享PED的计算结果l_jkd_k，所以可以将此结果缓存于先入先出缓存器(first in, first out, FIFO)中的直接用于l_ij的计算，从而避免重复运算，节省资源。PED的输入l_jk通过四选一数据选择器，按照一定的顺序将4组存储矩阵L中元素的随机存取存储器(random access memory, RAM)中的数据读出，用于计算d_j。矩阵A的数据通过数据分配器从外部传输到PEL与PED。其中，矩阵中a_ij与d_j, l_ij各自关联，但是PED与PEL计算存在延迟，故将数据存于FIFO_a₀~FIFO_a₄中。当前l_ij的计算结果与当前d_j值关联很大，且l_ij的计算还会用到计算d_j时的中间值l_jkd_k，所以PED要在PEL之前执行。在理想的情况下，为了实现运算速度最快，PEL的数量越多，效率越高。输入矩阵A为5×5的矩阵，在每一个内层循环中，需要计算的矩阵D的对角线元素d_j仅为1个，采用一个PED单元即可；而计算的矩阵L的元素l_ij的个数最多为4个，且互不依赖，故采用4个PEL并行计算。分解模块具体结构如图 3所示。其中RAMD为存储矩阵D中元素的随机存取存储器。

Figure 3.  Structure diagram of decomposition module

图 3中，PED和PEL模块均采取流水化设计方案，运算单元主要包括乘法器(multiplication)、除法器(division)、加法器(addition)和减法器(subtraction)，均使用Quartus Ⅱ提供的IEEE 754规范的浮点数运算IP核实现。

在设计PED和PEL模块时，鉴于乘法器在计算时存在固有延时，所以将输入的数据l_jk写入FIFO₁中实现缓存，此时FIFO₁就作为第2个乘法器的输入端，与第1个乘法器的计算结果进行第二次乘法运算。由于l_jkd_k为PED和PEL计算的共有部分，为减少乘法器的使用，故将计算中间结果l_jkd_k存于FIFO_ld中。累加器是由加法器和锁存器组成，每隔加法器计算的延时时间就将FIFO₂中数据读出一个，累加结束后，将结果传递给减法器输入端，并使能减法器。最终的计算结果存入RAM中。内层循环结束后，循环计算5次，即可获得对角阵D和下三角阵L矩阵的元素d_j和l_ij。

分解结束后，获得下三角矩阵L和对角阵D，根据(LDL^T)x=b，可以等价成3个简单的方程组Lz=b，Dr=z和L^Tx=r，由此可得求解公式:

式中，l_mp和l_um分别表示矩阵L中第m行p列和第u行m列所对应的元素；d_m表示矩阵D对角线上对应的元素；z_m, r_m, x_m, b_m分别代表列向量z, r, x, b的第m个元素；z_p，x_u分别代表列向量z的第p个元素和列向量x的第u个元素；q为待求向量的维数(q=5)。

求解模块要用到PEL模块中计算的数据l_ij，故每算出一个l_ij的值，求解模块中就要有一个计算单元与之相对应。求解模块的结构框图与PED和PEL的框图基本相同。

本文中基于Quartus Ⅱ 14.0的开发环境，将上述设计在Altera公司的DE1-SOC开发板上进行了验证。鉴于作者无法准确地得到光斑图像真实的中心，因此通过MATLAB程序生成一幅人工光斑图像，并附着一些噪声及干扰，该光斑图像的理论中心坐标为(150，166)，分别在FPGA平台和计算机Visual Studio(VS)环境下利用本文中的算法对光斑进行中心定位，根据(3)式计算的最终结果对比如表 1所示。其中，FPGA的工作时钟设定为50MHz，并采取IEEE 754规范下的单精度浮点制IP核来实现设计中的各浮点运算。

由表 1可知，本方法对圆的中心定位的精度小于0.1pixel，但FPGA平台与计算机VS环境下实现相同的算法计算的结果有差异，出现此结果的原因是，在FPGA平台使用单精度浮点IP核进行计算时，首先要将数据转换为单精度浮点数(32bit)，但对于较大的数据，尾数部分会有部分省略，进而对结果产生一定的误差。不过两个平台计算的中心位置的误差都小于0.1pixel，所以该方法保证了中心定位算法的精度。

保证中心定位算法精度的前提下，为测试方法的效率，将实现LDL^T解最小二乘问题的方法在FPGA平台和计算机VS下进行了计算效率比较。具体结果如表 2所示。

由表 2可知，在FPGA平台下实现解5维最小二乘问题时，本文中方法在保证精度的前提下，可达到相对于计算机21倍以上的速度。

光斑中心位置检测是光纤光束测量系统中的一个关键步骤。针对现有的基于计算机平台实现光斑参量测量实时性差，数据吞吐量不高的问题，本文中提出一种基于FPGA光纤光斑中心定位的方法：先对光斑图像进行中值滤波、二值化和边缘提取等预处理，提高算法的抗干扰能力，而后在FPGA平台上将获得的光斑边界数据进行分析处理。实验结果表明，在单精度浮点制条件下，本文中方法对光斑中心定位的绝对误差小于0.1pixel，在保证精度的同时可以实现21倍以上的速度提升。