CN103162818B - 基于矩不变性的激光光束束宽评价方法 - Google Patents

Abstract

本发明公开了一种基于矩不变性的激光光束束宽评价方法，属于光学计量领域。本评价方法首先对采集到的光束图像进行扣除背景、二值化、中值滤波处理得到滤波图像，之后在滤波图像上进行区域标记找出实测光斑区域，然后计算实测光斑的质心和二阶矩并将计算结果代入基于矩不变性原理获得的公式中，计算出光斑的长轴和短轴，即束宽；最后用与计算结果相适应的矩形框住实测光斑并显示。用本评价方法不仅可以计算规则形状，如圆形、椭圆形的激光光斑，也可以计算不规则形状的激光光斑，且计算结果不因光斑的平移和旋转而改变，反应了激光光束的固有特性。

Description

基于矩不变性的激光光束束宽评价方法

技术领域

本发明属于光学计量领域，主要涉及一种激光空域参数测量算法，尤其涉及一种激光光束束宽的评价方法。

背景技术

近年来，激光技术获得了突飞猛进的发展，在工业、医疗、国防、信息科学、生物技术等诸多领域都得到了广泛的应用。光束质量是判断激光系统综合性能的一项重要技术指标，对激光光束质量进行准确测量可以为激光输出提供准确可靠的实验数据，从而使激光更安全、可靠、有效地应用于各个领域。随着激光技术的发展，激光光束质量测量技术也在不断发展之中。

激光光束质量主要包括激光光束束宽、空间强度分布、束散角和指向稳定性等空域参数，这些参数决定着激光的有效传输距离、激光在目标上的聚焦性能等指标。在反映激光光束质量的空域参数中，激光光束束宽具有重要的地位，并直接决定着束散角的测量结果。在实际测量中，一般是通过阵列探测器(如CCD)对激光光束的横截面进行成像，再对采集的二维图像进行处理得到激光光束的束宽。

针对束宽的算法主要有1/e²定义法和二阶矩定义法。

1/e²定义法是经常被使用的一种方法，对旋转对称光束，在通过光斑质心的光强(光功率度)分布曲线I(r)上，最大值_max的1/e²(见如下公式1)处二点间距离之半定义为束宽w_e：

I(w_e)＝1/e²×I_max        (1) 

二阶矩定义法主要基于1995年ISO/D IS11146国际标准文件规定的束宽 定义及其测定方法。即用功率密度分布的二阶矩给出x，y方向的束宽w_x和w_y(见如下公式2和3)，上述积分原则上是对整个x-y平面积分，实际上此积分是至少在捕获光功率(能量)99％以上的区域进行。此方法在理论上较为严格，是当前广为应用的计算方法。

${w_y}^2=\frac{4{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}{(y-\stackrel{\‾}{y})}^{2}I(x,y,z)\mathrm{dxdy}}{{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}I(x,y,z)\mathrm{dxdy}}---\left(3\right)$

在实际的应用中，激光光斑不是严格服从高斯分布的旋转对称图形，如激光引信中应用的激光光斑由三个激光光束叠加而成，其形状类似于长条形的豆荚，还有的激光光斑呈现一些不规则的形状，如图1所示，这时采用1/e²定义法无法计算，而采用二阶矩定义法存在局限性，此局限性表现在，对具有一定光强分布的非旋转对称的光斑，其定义的x，y方向的束宽会随着光斑的旋转发生变化，即对同一光斑，当光斑在平面中发生旋转时，计算出的束宽会发生较大的变化，图2和图3直观显示了这种变化。图2是一束椭圆形激光光斑，其倾角为0°，即其长轴和短轴与坐标系的横轴和纵轴平行，采用二阶矩定义计算的束宽与椭圆形光斑的长轴和短轴相等；图3中同一椭圆光斑在CCD像面中旋转了45°，其长轴与CCD像面的横轴的夹角为45°，此时采用二阶矩定义计算的束宽与刚才计算出的0°倾角光斑的束宽存在倍的偏差。

发明内容

本发明要解决的技术问题是，针对现有技术的不足，提供一种以矩不变性原理为核心的激光光束束宽的评价方法。

为解决上述技术问题，本发明提供的评价方法包括以下步骤：

第一步，接到键盘命令后读取光束图像O和背景噪声图B，将光束图像O的灰度值减去背景噪声图像B中对应像素的灰度值，得到待测图像G并存储到计算机的硬盘中；

第二步，统计得到待测图像G的灰度直方图，计算所述灰度直方图的门限阈值m，根据下式处理得到所述待测图像G的二值化图像G₂，并将该二值化图像G₂存储在计算机的硬盘中：

$g_2(s,t)=\left\{\begin{array}{c}0\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;g(s,t)<m;\\ 255\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;g(s,t)\&GreaterEqual;m;\end{array}\right.$

其中，g₂(s，t)是二值化图像G₂在像素点(s，t)的灰度值；

第三步，用3×3像素模板对二值化图像G₂进行中值滤波处理，即把所述模板整体放置在二值化图像G₂的左上角，对模板所覆盖下的二值化图像G₂的9个灰度值进行大小排序，并将第五个值赋给中心像素，逐行逐像素移动模板，每移动一个像素就进行一次中值滤波，直至遍历二值化图像G₂的所有区域；

第四步，找出滤波图像G₃中灰度值为255的K个连通区域，并用q_k表示连通区域，其中，k＝1、2、3……K；

第五步，计算各连通区域所含的像素数，将具有最多像素数的连通区域q_k视为光斑区域；

第六步，将滤波图像G₃上的光斑区域对应到待测图像G上，获得光斑区域所对应的实测光斑Q；

第七步，根据以下步骤计算实测光斑Q的相关参数：

7.1根据下式计算实测光斑的质心坐标

$\begin{array}{c}\stackrel{\&OverBar;}{x}=\frac{\mathrm{\&Sigma;\&Sigma;s}\&times;g(s,t)}{\mathrm{\&Sigma;\&Sigma;g}(s,t)}\\ \stackrel{\&OverBar;}{y}=\frac{\mathrm{\&Sigma;\&Sigma;t}\&times;g(s,t)}{\mathrm{\&Sigma;\&Sigma;g}(s,t)}\end{array}$

其中，s是实测光斑中任一像素点的横坐标值，t是实测光斑中任一像素点的纵坐标值，(s，t)的取值范围为实测光斑Q的所有区域；

7.2通过下式计算实测光斑Q的长半径a、短半径b以及长半径与待测图像G水平方向的倾角θ：

$\left\{\begin{array}{c}\mathrm{\&theta;}=(1/2){\tan }^{-1}[2u_{11}/(u_{20}-u_{02})]\\ a={\left[\frac{u_{20}+u_{02}+{[{(u_{20}-u_{02})}^2+4{u_{11}}^2]}^{1/2}}{u_{00}/2}\right]}^{1/2}\\ b={\left[\frac{u_{20}+u_{02}-{[{(u_{20}-u_{02})}^2+4{u_{11}}^2]}^{1/2}}{u_{00}/2}\right]}^{1/2}\end{array}\right.$

$u_{\mathrm{pq}}=\&Integral;\&Integral;{(s-\stackrel{\&OverBar;}{x})}^p{(t-\stackrel{\&OverBar;}{y})}^{q}g(s,t)\mathrm{dsdt}$

其中，u_pq为实测光斑Q的二阶中心矩且p＝0、1、2，q＝0、1、2，积分范围为实测光斑Q包含的所有像素区域；

第八步，计算实测光斑Q的束宽，即长轴L_长＝2a，短轴L_短＝2b；

第九步，在屏幕上显示待测图像G、激光光斑的长轴L_长、短轴L_短以及长轴与待测图像G水平方向的倾角θ；并用长边等于长轴L_长，短边等于短轴L_短的矩形框套住所述实测光斑Q，其中，所述矩形框的中心与实测光斑Q的质心重合，矩形框的长边与待测图像G水平方向的夹角等于所述倾角θ。

本发明的技术效果体现在：

1.本发明在通过一系列搜索算法找到光斑区域后，计算光斑的质心和几个二阶矩，并进一步利用此质心和二阶矩的计算结果计算光斑的长短轴，计算过程不依赖光斑的形状，因此不仅可以计算光斑为圆形和椭圆形的激光光束，还可以用来计算其它不规则形状的光斑，并可以给出光斑长轴、短轴与坐标轴的夹角；

2.本发明的束宽评价方法基于矩的不变性原理设计，公式推理和计算结果均表明，对具有一定空间强度分布的光斑，当其在像面范围内发生平移和旋 转时，计算得到的光束束宽不变，因此本发明的评价结果反应了光束的固有特性参数。

附图说明

图1是不规则形状激光光斑的示意图。

图2是采用二阶矩定义法计算0°倾角的椭圆光斑束宽示意图。

图3是采用二阶矩定义法计算45°倾角的椭圆光斑束宽示意图。

图4是本发明评价方法的实现流程图。

图5是采用本发明评价方法计算45°倾角的椭圆光斑束宽示意图。

图6是采用本发明评价方法计算不规则形状激光光斑束宽示意图。

具体实施方式

下面结合附图及优选实例对本发明作进一步的详述。

本发明通过一台装有光束质量评价模块的计算机完成，前提条件是计算机硬盘中存储有光束横截面所成的光斑图像O(以下简称光束图像O)和背景噪声图像B，二者的分辨率均为M×N，即行数为M，列数为N。当计算机开机后，光束质量评价模块按照图4所示的流程执行以下操作步骤。

第一步，扣除背景，获得待测图像G

接到键盘命令后读取光束图像O，得到表示激光光束光强分布的二维灰度矩阵；再读取背景噪声图B，将光束图像扣除背景噪声，即将光束图像O的灰度值减去背景噪声图像B中对应像素的灰度值，得到待测图像G并存储到计算机的硬盘中，其具体算法公式如下：

g(s，t)＝o(s，t)-b(s，t)且s＝1、2、3......M，t＝1、2、3......N

式中，o(s，t)、b(s，t)依次为光束图像O、背景噪声图像B在像素点(s，t)位置的灰度值，g(s，t)为待测图像G在像素点(s，t)位置的灰度值。当整帧图像预处理 完毕后，就获得了待测图像G的二维灰度矩阵。

第二步处理待测图像G得到对应的二值化图像G₂

统计得到待测图像G的灰度直方图，并利用通用的“最小错误分割法”计算获得待测图像G灰度直方图的门限阈值m。根据下式处理得到二值化图像G₂，并将二值化图像G₂存储在计算机的硬盘中：

$g_2(s,t)=\left\{\begin{array}{c}0\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;g(s,t)<m;\\ 255\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;g(s,t)\&GreaterEqual;m;\end{array}\right.$

其中，g₂(s，t)是二值化图像G₂在像素点(s，t)的灰度值。

第三步，对二值化图像G₂进行中值滤波处理

实际采集的图像往往存在噪声，经过二值化处理后的噪声影响类似离子斑闪烁，如果不进行中值滤波，在下一步进行区域标记时会增加计算量，因而需要在此步中剔除离散的噪声。

中值滤波的过程是首先将如下3×3像素阵列放在二值化图像Gx的左上角，此阵列刚好覆盖二值化图像G₂的9个像素；

g2(s，t)	g2(s，t+1)	g2(s，t+2)
			g2(s+1，t)	g2(s+1，t+1)	g2(s，t+2)
g2(s+2，t)	g2(s+2，t+1)	g2(s+2，t+2)

其中，s的变化范围是1～M-2，t的变化范围是1～N-2；

对模板所覆盖下的二值化图像G₂的9个灰度值进行大小排序，并将第五个值(即中值)赋给中心像素，即g2(s+1，t+1)，边缘像素保留原来的值，之后所述3×3像素阵列在二值化图像G₂上逐行逐像素移动直至遍历二值化图像G₂的所有区域，每次移动都重复所述排序-赋值的中值滤波处理。处理完毕后，得到滤波图像G₃。

在此步骤中图像最边缘的像素没有处理，这不会对本发明产生影响，因为光斑成像的区域不会覆盖边缘。

第四步，在滤波图像G₃上进行区域标记

找出滤波图像G₃中灰度值为255的K个连通区域，并将全部连通区域标记为q_k，其中，k＝1、2、3......K。

第五步，找出滤波图像G₃中的二值化光斑

计算各连通区域所含的像素数，因为光斑对应的连通区域远远大于由噪声等干扰因素对应的连通区域，故将具有最多像素数的连通区域q_k视为光斑区域。

第六步，得到实测光斑区域Q

将滤波图像G₃上的光斑区域对应到待测图像G上，获得光斑区域所对应的实测光斑Q。

第七步，计算实测光斑Q的相关参数 

7.1根据下式计算实测光斑的质心坐标

$\stackrel{\&OverBar;}{x}=\frac{\mathrm{\&Sigma;\&Sigma;s}\&times;g(s,t)}{\mathrm{\&Sigma;\&Sigma;g}(s,t)}$

$\stackrel{\&OverBar;}{y}=\frac{\mathrm{\&Sigma;\&Sigma;t}\&times;g(s,t)}{\mathrm{\&Sigma;\&Sigma;g}(s,t)}$

其中，s是实测光斑中任一像素点的横坐标值，t是实测光斑中任一像素点的纵坐标值，(s，t)的取值范围为实测光斑Q的所有区域。

7.2计算实测光斑Q的长半径、短半径和倾角θ

通过下式计算实测光斑Q的长半径a、短半径b以及长轴与图像水平方向的倾角θ(以下简称倾角θ)。

$\left\{\begin{array}{c}\mathrm{\&theta;}=(1/2){\tan }^{-1}[2u_{11}/(u_{20}-u_{02})]\\ a={\left[\frac{u_{20}+u_{02}+{[{(u_{20}-u_{02})}^2+4{u_{11}}^2]}^{1/2}}{u_{00}/2}\right]}^{1/2}\\ b={\left[\frac{u_{20}+u_{02}-{[{(u_{20}-u_{02})}^2+4{u_{11}}^2]}^{1/2}}{u_{00}/2}\right]}^{1/2}\end{array}\right.$

$u_{\mathrm{pq}}=\&Integral;\&Integral;{(s-\stackrel{\&OverBar;}{x})}^p{(t-\stackrel{\&OverBar;}{y})}^{q}g(s,t)\mathrm{dsdt}$

其中，u_pq为实测光斑Q的二阶中心矩且p＝0、1、2，q＝0、1、2，积分范围为实测光斑Q包含的所有像素区域。当p＝0，q＝0时，即u₀₀的物理含义是实测光斑Q的面积，它与实测光斑Q在图像中的位置无关。

矩是一种统计特征，在图像处理中常用不同阶次的矩计算对象的位置、方向和尺度变换参数。根据矩的不变性原理，以下两个参数I₁和I₂具有在图像平移、旋转、缩放时计算结果保持不变的性质。对于本发明而言，I₁＝u₂₀+u₀₂，将I₁和I₂带入a、b的计算式中，可得：

$\left\{\begin{array}{c}a={\left[\frac{I_1+I_2^{1/2}}{u_{00}/2}\right]}^{1/2}\\ b={\left[\frac{I_1-I_2^{1/2}}{u_{00}/2}\right]}^{1/2}\end{array}\right.$

上式表明：实测光斑Q的长半径a、短半径b的计算结果既不会因为实测光斑Q在图像中的平移而改变，也不会因为实测光斑Q的旋转而改变，即具有旋转不变性，因此反应了实测光斑的固有参数。

第八步计算实测光斑Q的束宽

计算实测光斑Q的长轴L_长＝2a，短轴L_短＝2b，并将L_长和L_短作为实测光斑的束宽。

对于形状规则的光斑，本发明计算结果同其几何参数完全一致。例如， 如果实测光斑为圆形，则长轴L_长和短轴L_短相等，并等于实测光斑的直径；如果实测光斑为椭圆，则长轴L_长与其长轴相等，短轴L_短与其短轴相等；对于其它如图1所示的不规则形状的光斑，同样可以将本方法的计算结果作为其束宽，并将倾角θ作为所述光斑与横坐标的夹角。

第九步，显示计算数值并图形化显示计算结果

在屏幕上显示待测图像G、激光光斑的长轴L_长、短轴L_短即光束束宽以及倾角θ的计算结果；并用长边等于长轴L_长，短边等于短轴L_短的矩形框住实测光斑区域，其中，矩形的中心与实测光斑Q的质心重合，矩形的长边与待测图像G水平方向的夹角等于倾角θ(参见图5、图6)。

Claims (1)

1.一种基于矩不变性的激光光束束宽评价方法，其特征在于：该方法包括以下操作步骤：

第一步，接到键盘命令后读取光束图像O和背景噪声图B，将光束图像O的灰度值减去背景噪声图像B中对应像素的灰度值，得到待测图像G并存储到计算机的硬盘中；

第二步，统计得到待测图像G的灰度直方图，计算所述灰度直方图的门限阈值m，根据下式处理得到所述待测图像G的二值化图像G₂，并将该二值化图像G₂存储在计算机的硬盘中：

$g_2(s,t)=\left\{\begin{array}{c}0......g(s,t)<m;\\ 255......g(s,t)\&GreaterEqual;m;\end{array}\right.$

其中，g₂(s，t)是二值化图像G₂在像素点(s，t)的灰度值；

第三步，用3×3像素模板对二值化图像G₂进行中值滤波处理，即把所述模板整体放置在二值化图像G₂的左上角，对模板所覆盖下的二值化图像G₂的9个灰度值进行大小排序，并将第五个值赋给中心像素，逐行逐像素移动模板，每移动一个像素就进行一次中值滤波，直至遍历二值化图像G₂的所有区域；

第四步，找出滤波图像G₃中灰度值为255的K个连通区域，并用q_k表示连通区域，其中，k＝1、2、3......K；

第五步，计算各连通区域所含的像素数，将具有最多像素数的连通区域q_k视为光斑区域；

第六步，将滤波图像G₃上的光斑区域对应到待测图像G上，获得光斑区域所对应的实测光斑Q；

第七步，根据以下步骤计算实测光斑Q的相关参数：

7.1根据下式计算实测光斑的质心坐标

$\stackrel{\&OverBar;}{x}=\frac{\mathrm{\&Sigma;\&Sigma;s}\&times;g(s,t)}{\mathrm{\&Sigma;\&Sigma;g}(s,t)}$

$\stackrel{\&OverBar;}{y}=\frac{\mathrm{\&Sigma;\&Sigma;t}\&times;g(s,t)}{\mathrm{\&Sigma;\&Sigma;g}(s,t)}$

其中，s是实测光斑中任一像素点的横坐标值，t是实测光斑中任一像素点的纵坐标值，(s，t)的取值范围为实测光斑Q的所有区域；

7.2通过下式计算实测光斑Q的长半径a、短半径b以及长半径与待测图像G水平方向的倾角θ：

$\left\{\begin{array}{c}\mathrm{\&theta;}=(1/2){\tan }^{-1}[{2u}_{11}/(u_{20}-u_{02})]\\ a={\left[\frac{u_{20}+u_{02}+{[{(u_{20}-u_{02})}^2+{{4u}_{11}}^2]}^{1/2}}{u_{00}/2}\right]}^{1/2}\\ b={\left[\frac{u_{20}+u_{02}-{[{(u_{20}-u_{02})}^2+{{4u}_{11}}^2]}^{1/2}}{u_{00}/2}\right]}^{1/2}\end{array}\right.$

$u_{\mathrm{pq}}=\&Integral;\&Integral;{(s-\stackrel{\&OverBar;}{x})}^p{(t-\stackrel{\&OverBar;}{y})}^{q}g(s,t)\mathrm{dsdt}$

其中，u_pq为实测光斑Q的二阶中心矩且p＝0、1、2，q＝0、1、2，积分范围为实测光斑Q包含的所有像素区域；

第八步，计算实测光斑Q的束宽，即长轴L_长＝2a，短轴L_短＝2b；

第九步，在屏幕上显示待测图像G、激光光斑的长轴L_长、短轴L_短以及长轴与待测图像G水平方向的倾角θ；并用长边等于长轴L_长，短边等于短轴L_短的矩形框套住所述实测光斑Q，其中，所述矩形框的中心与实测光斑Q的质心重合，矩形框的长边与待测图像G水平方向的夹角等于所述倾角θ。