CN109859135A - 一种应用于关联成像的图像增强处理方法 - Google Patents

Abstract

本发明公开了一种应用于关联成像的图像增强处理方法，包括：对参考臂光场采样数据与信号臂光强响应数据进行关联计算，得到关联成像结果图；计算结果图归一化自相关矩阵、半高全宽；将结果图减去归一化自相关矩阵和增益系数的乘积得到残差图；迭代计算残差图内最大值点及对应的归一化自相关矩阵，并用残差图中减掉该矩阵与增益系数的乘积；将该点的位置与强度计入净点图中；满足迭代终止条件时停止迭代；用半高全宽拟合高斯函数与净点图卷积，得到无周期性结构的关联成像结果图。本发明可解决空间频率采样不充分造成的周期性重复结构问题，提升关联成像结果图的可见度和信噪比，该方法抗噪声能力强，计算速度快，实用效果好。

Description

一种应用于关联成像的图像增强处理方法

技术领域

本发明属于图像处理领域，特别涉及一种基于关联成像中二阶自相关函数的关联成像结果图像去周期性重复结构的图像增强处理方法。

背景技术

关联成像技术是基于照明电磁场特定阶数的相关性质进行的一种非直接多次曝光成像方式。其对物探测的光路无需空间分辨的“离物成像”特性大大降低了成像系统对面阵探测器的要求，能够实现特定条件下的超衍射极限空间分辨，并能实现与激光雷达、三维立体成像等应用光学和成像领域技术的高效结合，具有重要的应用意义与广泛的应用前景。

目前关联成像技术的主流应用方向之一是结合激光雷达的技术特点，实现长距离、高空间分辨性能的关联成像激光雷达。通常情况下，关联成像激光雷达因功率、造价等因素制约，发射端的光纤激光头数量有限，对应空间频率采样不足以覆盖整个空间频率域，导致成像系统的点扩散函数有周期性延拓的旁瓣，成像结果中产生周期性重复结构，严重影响图像内容的可见度。因此，需要在保留图像中的有效信息的同时，采用图像处理方法将图像中的空间周期性重复结构去除，提升图像内容的信噪比以及可见度。

目前可用于关联成像的图像去周期性重复结构的方法有压缩感知算法、矩阵伪逆算法等。压缩感知算法由E.Candès，D.Donoho和T. Tao提出，基于图像的稀疏先验特性实现突破奈奎斯特采样极限的图像重构恢复算法。应用该算法需要提供成像目标在某个域上的稀疏先验特性，而往往在对目标进行成像之前该特性并不能确知。而关联成像的矩阵伪逆算法由C.Zhang等人提出，基于对重复测量矩阵的伪逆计算实现更高信噪比、更高分辨率的关联图像恢复。该方法的抗噪声性能差，矩阵伪逆计算的时间复杂度高，成像速率难以适应关联成像激光雷达的实际应用场景。

与关联成像激光雷达所面临的空间频率采样不充分问题类似，现代天文观测中所使用的甚大望远镜阵列、甚长基线干涉阵列等技术的空间频率采样同样不充分，天体的点源响应也有彼此覆盖的旁瓣，影响天体图像的可见度。而天文图像处理领域已发展出一套成熟的 CLEAN数值反卷积算法增强天体图像的可见度。但与关联成像激光雷达所不同的是：其一，关联成像激光雷达系统中此前并未有对应于天文成像系统点源响应的概念，CLEAN算法很难直接应用于关联成像激光雷达系统；其二，CLEAN算法需要计算每幅残差图内最大值点处的点源响应，算法的时间复杂度高，同时对内存空间有较高的消耗。

发明内容

本发明是为了提供一种基于关联成像中的高阶自相关函数的去周期性结构的图像增强处理方法，该算法能够解决关联成像结果图中由于空间频率采样不足产生的图像周期性结构问题，提高了成像结果的恢复质量并降低了图像背景噪声。

相较于其他的关联成像去图像周期性结构的方法无需探测目标的稀疏先验特性约束，抗噪声性能强，时间复杂度低，同时极大地降低了算法对内存空间的消耗。相较于传统CLEAN算法，改进后的本方法无需在每一次迭代中都遍历全部的光场信息并存储迭代结果，仅需在第一次迭代运算中遍历全部的光场信息并计算一次中间结果即可，时间复杂度从o(n²)降至o(n)，同时算法运行占用的内存空间也显著降低。

本发明提供的技术方案为：

一种基于关联成像中的高阶自相关函数的去周期性结构的图像增强处理方法，包括：

关联成像的参考光路光场空间分布和信号光路光强探测响应两类重复测量的数据，其中参考光路光场分布的重复测量数据以视频的形式存储，信号光路光强响应以文本形式存储。

首先，进行关联成像矩阵预处理：

将参考臂光场采样数据与对应的信号臂光强响应数据进行关联运算，获取初始的关联成像结果图，将其作为本增强方法的输入，记为脏图；

初始化残差图，残差图初始值等于脏图；

找到脏图中强度最大值的点的对应位置，并记录该点的坐标与强度值；

计算每幅参考臂光采样数据中对应脏图最大值坐标的点的值与当前参考臂光场分布矩阵的乘积，并按采样数进行累加，记为该点处的自相关矩阵；

将每幅参考臂采样数据中对应脏图最大值坐标点的值，存为最大值数组；

将每幅参考臂采样数据分布累加，存为参考臂累加矩阵；

自相关矩阵除以参考臂累加矩阵与最大值数组均值的乘积，记为归一化自相关矩阵；

创建新矩阵g2,其大小等于归一化自相关矩阵的四倍；

将归一化自相关矩阵平移，使归一化自相关矩阵最大值点与g2中心点重合后，将归一化自相关矩阵赋值给g2矩阵对应位置；

将g2矩阵中其他未赋值位置统一取值为归一化自相关矩阵的像素平均值；

下一步，开始进行迭代去周期结构算法：

初始化全零的净点图矩阵，矩阵大小等于归一化自相关矩阵；

设定迭代次数上限、迭代增益和迭代终止比值，其中迭代终止比值为残差图最大值比值和脏图最大值比值，开始循环：

找到当前残差图中最大值点，计算残差图中最大值点坐标与脏图中最大值点坐标之间的偏移量；

在净点图矩阵中，将残差图中最大值点的对应位置赋值为残差图最大值与迭代增益的乘积；

以残差图中最大值点坐标为中心，在g2矩阵中截取与归一化自相关矩阵相同大小的矩阵，记为脏束；

将残差图减掉脏束与迭代增益的乘积以更新残差图；

判断残差图中像素的最大值是否小于迭代终止比值；

若是残差图中像素的最大值小于迭代终止比值或循环次数大于迭代次数上限，跳出循环，否则继续循环；

选出g2矩阵最大值点处的一行像素值，对其根据实际情况进行平移处理后，另存为g2Lane数组；

找到g2Lane数组的最大值，将该数值除以2作为g2Lane的半高，并将整个g2Lane数组减去该半高值，并取结果的绝对值作为新的g2Lane数组；

取与x轴重合的水平直线作为控制变量，设定每次上移该控制变量的步长值，并按该步长值向上平移控制变量，并统计控制变量与新的g2Lane数组相交点的个数与其横坐标值；

当交点的个数恰好等于二且横坐标值的差的绝对值大于一时，结束半高全宽计算过程，否则向上平移控制变量；

归一化自相关函数的半高全宽等于两个交点横坐标之差的绝对值。

拟合二维高斯函数矩阵，矩阵大小等于归一化自相关矩阵，高斯函数的标准差等于半高全宽除以一个特定值，该值可根据经验确定；

下一步，进行净点集合去噪与图像重构：

设定邻域半径与最小点数阈值，设定全零的矩阵作为物点集合矩阵，其大小等于净点图矩阵；

记录净点图中所有非零值的对应坐标，保存为数组净点序列；

遍历访问净点序列，找到净点序列当前值在净点图矩阵中的对应点坐标；

统计以该点坐标为中心，欧氏距离小于邻域半径的非零点个数；

若非零点个数小于最小点数，该点判定为噪声点；否则判断为物点，将物点坐标在物点集合矩阵中的对应处赋值为一；

遍历完成后，将净点图矩阵与物点集合矩阵相乘，更新净点图并与所拟合的二维高斯函数矩阵进行二维卷积，存为净图；

净图即为无周期性结构的关联成像结果图。

本发明至少包括以下有益效果：针对关联成像图像中具有空间周期性结构的问题，首先应用了CLEAN算法进行图像周期性结构去除，并针对于算法的不足之处，从算法与关联成像系统结合与时间复杂度的角度提出了改进，一是在算法中利用关联成像参考光路的自相关矩阵替代原CLEAN算法中的脏束，保证了CLEAN算法在关联成像系统中的有效应用；二是改变了原CLEAN算法中每个不同点对应的脏束均要重新计算的低效方法，利用初次脏图最高点对应的自相关矩阵赋值入双倍大小的g2矩阵，并以不同残差图内最大值点为中心的平移提取矩阵加以替代，显著降低了算法的计算量和运行时间，极大的提升算法的实用性，使得算法足以完成图像的实时重构工作；附加了利用密度聚类方法对图像恢复结果的去噪处理方案，显著减少了原 CLEAN算法常见的噪点数量，提升了恢复图像的信噪比以及可见度。

附图说明

图1是实施例一的图像增强处理方法(部分)的流程图；

图2是实施例二的图像增强处理方法的流程图；

图3是g2矩阵计算流程图；

图4是实施例一的预期目标物体图；

图5是实施例一的二阶关联结果图；

图6是实施例一的图像增强处理方法(部分)的结果图；

图7为实施例二中预期目标物体图；

图8为实施例二中二阶关联运算的结果图；

图9为实施例二中传统CLEAN算法的结果图；

图10为实施例二中图像增强处理方法的结果图；

具体实施

下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

实施例一与实施例二的数据均来源于合作单位的外场试验，在获得所采集数据后对数据进行关联计算后得到关联结果，在实例一中，只进行了本专利算法的部分步骤，以测试本专利算法对空间结构噪声的消除效果，在实例二中则进行了本专利算法的完整步骤，在完成去噪声、降低了算法运行所需时间的同时提高了图像优化质量。参考附图说明，可以看出，在原始的二阶关联运算结果图中难以分辨出实际物体的形状，背景可见周期性干扰噪声。通过实施例一可以看出本算法在去周期性结构噪声方面具有优秀的效果，通过实施例二可以看出本算法运行完整步骤后具有更加优秀的成像效果、更低的算法运行时间。

实施例一对应了权利要求2，3，5，7；实施例二对应了权利要求 2～7。

实施例一：

本实施例中参考臂光场采样数据包含有10000张图片，每张图片有224*224个像素，存储于avi格式中，另有信号臂光强响应信息，为一1*10000的向量数据，存储于文本文件中。在应用图像增强处理方法算法时，首先借助附加库读取视频文件，并将其内容存储于计算机内存中，之后获取信号臂光强响应信息，同样存储于内存之中。

附图说明中图1为本实施例的流程图，

附图说明中图3为本实施例中g2矩阵的计算流程图，

如图1所示，本实施例提供的应用于关联成像的图像增强化处理方法(部分)，包括如下步骤：

设定本专利算法的循环终止条件为：

当前残差图矩阵最大值为初始脏图矩阵最大值的1/4时跳出循环，或者循环次数达到10000次时跳出循环；

设定图像增强方法循环迭代增益为0.005；

初始化所需的各类矩阵为零矩阵；

进行关联计算：

在一次二阶关联的循环中：

依次获取视频文件中的参考臂光场采样数据并修改数据类型，保存为光场的矩阵数据；

累加视频帧数据得到全部参考臂光场采样数据的累加矩阵；

通过每一帧光场的矩阵数据与对应帧信号臂光强响应数据的运算得到关联矩阵；

计算归一化二阶强度关联；

结束该次二阶关联循环；

通过关联运算得到关联结果；

进行图像增强算法：

取二阶关联结果为初始残差图和脏图；

取每一帧光场矩阵算其最高值，并记录最高值坐标并计算归一化自相关矩阵；判断是否满足终止条件，若满足则跳出循环，反之进行后续步骤；

如图3流程图所示计算g2矩阵；

残差图减去当前帧最大值点对应的归一化自相关矩阵与增益系数的乘积，在净点集合中对应该最大值点的位置加入该最大值与增益系数的乘积以更新残差图，将当前帧当前最高值点加入净点图以更新净点图；

更新循环数，判断是否满足终止条件，若满足则跳出循环，反之进行后续步骤；

结束一次循环并进行下一次循环；

将净点图中每一点与计算得到的点扩散函数进行卷积运算，得到最终的无周期性结构的结果图。

附图说明中的图4为此实施例中目标物体预期的结果图，图5为此实施例中的二阶关联运算结果，其可见度为0.0021，信噪比为 0.1582；图6为此实施例中的增强结果，并计算了其可见度为0.9766，信噪比为3.3722，见表二。

实施例二：

本实施例中参考臂光场采样数据包含有10000张图片，每张图片有250*250个像素，存储于avi格式中，另有信号臂光强响应数据，为一1*10000的向量数据，存储于文本文件中。在应用本专利算法时，首先借助附加库读取视频文件，并将其内容存储于计算机内存中，之后获取光强信息，同样存储于内存之中。

附图说明中图2为本实施例的流程图，

附图说明中图3为本实施例中g2矩阵的计算流程图，

如图1所示，本实施例提供的应用于关联成像的图像增强处理方法(部分)，包括如下步骤：

设定图像增强化方法的终止条件为：

当前残差图矩阵最大值为初始脏图矩阵最大值的1/5时跳出循环或者循环次数达到10000次时跳出循环；

设定图像增强方法循环迭代增益为0.005；

设定聚类半径；

初始化所需的各类矩阵为零矩阵；

进行关联计算：

在一次二阶关联的循环中：

依次获取视频文件中的参考臂光场采样数据并修改数据类型，保存为光场的矩阵数据；

累加视频帧数据得到全部参考臂光场采样数据的累加矩阵；

通过每一帧光场的矩阵数据与对应帧信号臂响应数据的运算得到关联矩阵；

计算归一化二阶强度关联；

结束该次二阶关联循环；

通过关联运算得到关联结果图；

进行图像增强方法：

取第一帧参考臂光场采样数据矩阵为初始残差图和脏图，计算其最高值，并记录最高值坐标，之后用Max来表述该点；

进行一次对参考臂光场采样数据的遍历循环，并记录在Max点处对应的值，并计算归一化自相关矩阵；

如图3流程图所示计算g2矩阵：

初始化500*500的零矩阵为g2，并将其中间部分赋值为计算所得的归一化自相关矩阵，并将其他位置的值赋值为自相关矩阵的像素平均值；

寻找g2矩阵最高值点，并抽取该点所在列行程一维向量并计算其半高全宽；

开始一次图像增强方法循环：

在每一帧中寻找当前帧最高值点；

判断是否满足终止条件，若满足则跳出循环，反之进行后续步骤；

计算该点与点Max的欧氏距离，若欧氏距离小于聚类半径则加入净点图，反之抛弃该点；

无论欧氏距离如何，在残差图迭代减去最大值点对应的归一化自相关矩阵与增益系数的乘积，在净点集合中对应该最大值点的位置加入该最大值与增益系数的乘积以更新残差图；

更新当前循环数，判断是否满足终止条件，若满足则跳出循环，反之进行后续步骤；

结束此次循环并进行下一次循环；

根据半高全宽拟合高斯函数并与净点图进行卷积，得到最终的结果图。

附图说明中的图7为此实施例中目标物体图，图8为此实施例中的二阶关联运算结果，其可见度为0.2405，信噪比为0.0388，图10 为此实施例中的增强结果，并计算了其可见度为0.5829，信噪比为 1.7046，见表二。

本专利算法可与传统CLEAN算法进行比较，传统CLEAN算法计算结果如附图说明中的图9，其可见度为0.4264，信噪比为0.0565，见表二。

表一实施例二耗时对比

计算方法	平均计算耗时
		二阶关联	8s
未优化的CLEAN算法	270s
		本专利算法	10s

表二实施例一、实施例二图像可见度以及信噪比

尽管本发明的实施方案已公开如上，但其并不仅仅，限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (7)

1.一种应用于关联成像的图像增强处理方法，其特征在于，包括：

对参考臂光场采样数据与信号臂光强响应数据进行关联计算，得到初始的关联成像结果图；

找到初始的关联成像结果图内最大值点，记录其坐标，计算该最大值点对应的归一化自相关矩阵，根据该矩阵计算归一化自相关函数的半高全宽并记录该值；

在关联成像图中迭代减去归一化自相关矩阵与增益系数的乘积，在净点集合中的该最大值点对应的位置加入该最大值与增益系数的乘积；

设定邻域半径和最小点数阈值用以对净点图进行过滤，统计净点图中有值的各点其邻域半径范围内净点的个数，若净点个数小于最小点数阈值，则在净点集合中去掉该点，否则保留；

利用归一化自相关函数的半高全宽拟合高斯函数，与过滤后的净点图进行卷积运算，得到最终的无周期性结构的关联成像结果图。

2.如权利要求1所述一种应用于关联成像的图像增强处理方法，其特征在于，得到关联成像结果的过程具体包括：

将探测器所记录的每一帧参考臂光场采样数据相加得到光场的叠加矩阵，将该矩阵按采样数计算平均，计算得到参考臂光场叠加矩阵的均值矩阵；

计算每一帧参考臂光场采样数据与该均值矩阵的差值，记为该帧采样对应的参考臂光场涨落矩阵；

计算信号臂光强响应数据的均值，记为光强均值，并将每一帧光强响应数据减掉光强均值，存为信号臂光强涨落向量；

设定涨落关联矩阵，并将全零的矩阵作为初始的涨落关联矩阵。计算每一帧参考臂光场涨落矩阵与对应的信号臂光强涨落向量对应值的乘积，将结果累加入涨落关联矩阵，累加完成后，将涨落关联矩阵除以参考臂光场采样均值矩阵与信号臂光强响应的光强均值乘积，得到初始的关联成像结果图。

3.如权利要求1所述一种应用于关联成像的图像增强处理方法,其特征在于，计算初始的关联成像结果图内最大值点与最大值点对应的归一化自相关矩阵，具体包括：

寻找初始的关联成像结果图中的强度的最大值，得到该最大值与该最大值点的坐标；

若找到的最大值点的数目不为一，则选择找到的第一个最大值用于后续的计算；

设定全零的矩阵为初始的自相关矩阵，将与初始的关联成像结果图最大值点坐标处对应的每一帧参考臂光场采样数据中强度与该帧参考臂光场采样数据矩阵相乘，并与自相关矩阵相加得到累加后的自相关矩阵；

将累加后的自相关矩阵除以每一帧对应初始的关联成像结果图中最大值坐标处的点的强度的均值与参考臂光场采样均值矩阵的乘积，得到归一化自相关矩阵。

4.如权利要求1所述一种应用于关联成像的图像增强处理方法，其特征在于，计算归一化自相关函数的半高全宽具体包括：

选出归一化自相关矩阵最大值点处的一行，另存这一行数据为g2Lane数组；

找到g2Lane数组的最大值，将该数值除以2作为g2Lane的半高，并将整个g2Lane数组减去该半高值，并取结果的绝对值作为新的g2Lane数组；

取与x轴重合的水平直线作为控制变量，设定每次上移该控制变量的步长值，并按该步长值向上平移控制变量，并统计控制变量与新的g2Lane数组相交点的个数与其横坐标值；

当交点的个数恰好等于二且横坐标值的差的绝对值大于一时，结束半高全宽计算过程，否则向上平移控制变量；

归一化自相关函数的半高全宽等于两个交点横坐标之差的绝对值。

5.如权利要求1所述一种应用于关联成像的图像增强处理方法,其特征在于，在初始的关联成像图增强处理过程中迭代减去每一帧对应于初始的关联成像图最大值处的点与归一化自相关矩阵与增益系数的乘积，在净点图中对应初始的关联成像图最大值处的点的位置加入该乘积，具体包括：

设定全零矩阵为g2矩阵,其大小等于归一化自相关矩阵的四倍；

将归一化自相关矩阵平移，使归一化自相关矩阵最大值点坐标为g2矩阵中心点坐标，将归一化自相关矩阵赋值给g2矩阵对应位置；

将g2矩阵中其他未赋值位置统一取值，为归一化自相关矩阵的平均值；

设定迭代次数上限、迭代增益和迭代终止比值，其中迭代终止比值为残差图最大值比值和脏图最大值比值，开始迭代循环：

找到残差图中最大值点，计算残差图中最大值点坐标与脏图中最大值点坐标之间的偏移量；

在净点图矩阵中，将残差图中最大值点的对应位置赋值为残差图最大值与迭代增益的乘积；

以残差图中最大值点坐标为中心，在g2矩阵中截取与归一化自相关矩阵相同大小的矩阵，记为脏束；

将残差图减掉脏束与迭代增益的乘积以更新残差图；

判断残差图中像素的最大值是否小于迭代终止比值；

若是残差图中像素的最大值小于迭代终止比值或循环次数大于迭代次数上限，跳出循环，否则继续循环。

6.如权利要求1所述应用于关联成像的图像增强处理方法，其特征在于，进行净点集合去噪与图像重构，包括如下步骤：

设定邻域半径与最小点数阈值，设定全零的矩阵作为物点集合矩阵，其大小等于净点图矩阵；

记录净点图中所有非零值的对应坐标，保存为数组净点序列；

遍历访问净点序列，找到净点序列当前值在净点图矩阵中的对应点坐标；

统计以该点坐标为中心，欧氏距离小于邻域半径的非零点个数；

若非零点个数小于最小点数，该点判定为噪声点；否则判断为物点，将物点坐标在物点集合矩阵中的对应处赋值为一；

遍历完成后，将净点图矩阵与物点集合矩阵相乘，更新净点图。

7.如权利要求1所述一种应用于关联成像的图像增强处理方法，其特征在于，利用归一化自相关函数的半高全宽拟合高斯函数，与净点图相卷积，得到无周期性结构的关联成像图，具体包括：

拟合二维高斯函数矩阵，矩阵大小等于归一化自相关矩阵，高斯函数的中心位于矩阵中心元素；

将二维高斯函数矩阵与净点图矩阵进行二维卷积运算，得到无周期性结构的关联成像结果图。