CN106952241B - 一种基于形态学方法及Meanshift算法的电磁图像分块方法 - Google Patents

Abstract

本发明是基于形态学方法及Meanshift算法的电磁图像分块方法，本发明首先对图像进行同态滤波抑噪，然后使用形态学方法对图像进行处理，从而使图像更易于进行Mean Shift聚类分割，然后计算图像粗糙度和像素值平均偏移作为聚类算法的阈值参数。接下了进行MS分割，并对分割后图像判断是否为同一信号区域，对于同一信号区域进行合并。对方法进行了仿真和实验验证，结果显示方法可以高效准确的对多区域多频率电磁图像进行图像抑噪并分割。

Description

一种基于形态学方法及Meanshift算法的电磁图像分块方法

技术领域

本发明涉及对多区域多频率电磁辐射源成像的噪声抑制及分区图像恢复算法，具体涉及电磁探测及图像处理领域。

背景技术

目前对电磁干扰源的检测手段存在很多缺点，如检测速度慢、难以全面检测等，而利用抛物反射面进行电磁成像的方法可以快速准确地检测干扰源的位置。电磁成像研究目前主要集中在被动毫米波成像领域，受天线尺寸、衍射受限等因素的影响，所成图像分辨率不高，相关的提高分辨率的超分辨算法成为了研究的热点。而对于电磁干扰源的检测成像，与毫米波成像相比波长更长，衍射受限更为严重，且所成图像中包含多种干扰源的点扩展函数，如何提高其分辨率使得干扰源能够清晰分辨成了新的难点。

发明内容

本发明技术解决问题：克服现有技术的不足，提出了一种基于形态学方法及Meanshift算法的电磁图像分块方法，可以使图像处理符合人眼对于亮度响应的非线性特性，避免了直接对图像进行傅立叶变换处理的失真。

本发明技术解决方案：在进行分区算法之前，先使用同态滤波对图像进行去噪处理，同态滤波是把频率过滤和灰度变换结合起来的一种图像处理方法，它依靠图像的照度/反射率模型作为频域处理的基础，利用压缩亮度范围和增强对比度来改善图像的质量。

接下来对去噪处理后的图像进行区域抽离。区域抽离的基础是进行图像的分割。简单的图像分割方法包括均匀分割与非均匀分割均。采用均匀分割的方法即对所成图像横向分割或纵向分割，规定所分区域的个数，将图像均匀的分割，这种方法显然是不适用的，无法准确找到各干扰源位置，并且易将一个干扰源分到不同区域。非均匀分割法则根据各干扰源位置和大小进行分割，可以准确找到分辨率相同的干扰源所在区域，但需要手动进行，耗时长且没有统一的分割方法。作为图像处理领域的热点问题，近年来已有许多有效的图像分割方法，而其中基于Mean Shift即均值漂移的方法适用于对电磁辐射源成像的分割。

但Mean Shift图像分区算法倾向于识别凸性分布，大小相近，密度相近的聚类，因此需要对原始电磁探测成像进行处理，使得待分割区域边缘尽量平滑，且为凸性分布，从而提高Mean Shift算法的精确度。经过分割后得到的图像有割裂同一区域或重叠的可能性，因此在Mean Shift算法计算之后增加对于重叠图像及错误割裂图像的判断，得到的最终分割图像能够较准确的分别出各个区域。

本发明包括以下步骤：

步骤1，获得未知辐射源的降质图像，对图像进行对数运算，将图像转换为照度和反射两个部分相加的形式，得到取对数后图像；

步骤2，对取对数后图像进行滤波处理，滤除图像的高频噪声，从而压低亮度范围，再对图像取指数，恢复电磁功率图像，得到去噪处理后图像；

步骤3，对去噪处理后图像进行sorbel边缘提取，并使用图像形态学方法对边缘提取图像进行膨胀腐蚀操作，从而将图像填充为信号区域为凸性边缘的整体二值图像；

步骤4，计算整体二值图像的粗糙度，粗糙度代表了图像中平均纹理的大小，作为MeanShift算法中距离聚类阈值hs，然后计算图像像素平均偏移值，此偏移值代表图像各像素点之间像素值变化情况，作为MeanShift算法中像素值聚类阈值h_r；

步骤5，根据步骤4中计算得到的距离聚类参数h_s和像素值聚类参数h_r对填充图像进行MeanShift算法进行图像分块；

步骤6，根据步骤5中的分块情况，将步骤2去噪后图像的不同的区域抽出信号区域，然后分块图像补充为原有图像大小，方便后面进行图像恢复；

步骤7，对分块得到的图像进行判断，如果是同一个信号分块区域则进行合并。

本发明的优点与积极效果在于：

(1)传统图像分块算法对于需要分块的每一块区域的边缘清晰度比较严格的要求，对于边缘比较模糊的图像的边缘识别以及图像分割效果比较差。而对于电磁成像得到的图像来说，由于点扩散函数的影响其电场强度的分布趋势接近于sinc函数，即在边缘处像素值变化缓慢，因此传统图像分块算法对于电磁成像得到的电场强度图像的分辨效果较差，分区识别的正确率很低。而先采用形态学方法处理后，电磁成像的信号区域可以变换为边缘比较明显的且基本为凸性边缘图像区域，更易之后MeanShift聚类算法进行计算。

(2)由于电磁成像系统为宽带电磁成像，随着电磁信号频率的降低，电磁成像系统点扩散函数主瓣宽度增加其在像平面上成像区域大小也逐渐增加，因此对于分块算法来说，其分块的聚类参数需要根据信号区域大小进行调整。在本发明中，计算图像的粗糙度和像素值平均偏移来作为聚类算法的阈值参数，这两个参数在衡量了图像的像素值区域大小，变化快慢等，在一定程度上可以体现图像的纹理情况，因此使用其作为聚类算法的参数在一定程度上可以达到宽带区域分割的目的。

附图说明

图1为本发明方法实现流程图；

图2为无噪声情况下不同频率多辐射源电磁图像；

图3为无噪声情况下不同频率多辐射源电磁图像进行sobel边缘检测图像；

图4为图像形态学方法处理后的二值图像

图5经过Meanshift分割算法后得到的单信号图像，其中a为信号1提取后图像，b为信号2提取图像，c为信号3提取图像；

图6电磁成像系统实验得到含噪情况下不同频率多信号源电磁图像；

图7经过同态滤波后的低噪图像；

图8经过Meanshift分割算法后得到的单信号图像，其中a为信号1提取后图像，b为信号2提取图像；

图9单信号图像经过LR超分辨处理，合并多个信号超分辨后图像得到最终高分辨率电磁图像。

具体实施方式

下面将结合附图和实施例对本方法作进一步的详细说明。

如图1所示，本发明针对于同一电磁成像图像存在不同频率辐射源的问题，提出一种使用形态学和MeanShift聚类算法相结合的宽带电磁图像分割算法。对取对数后图像进行滤波处理，滤除图像的高频噪声，从而压低亮度范围，再对图像取指数，恢复电磁功率图像。

步骤1，获得未知辐射源的降质图像，对图像进行对数运算，将图像转换为照度和反射两个部分相加的形式。

即图像f(x,y)可以表示为照度i(x,y)和反射r(x,y)两部分的乘积：

f(x,y)＝i(x,y)r(x,y)

因此对于原始图像不能分别对照度和反射部分分别进行操作，因此选择对图像进行对数运算，假设z(x,y)＝lnf(x,y)＝lni(x,y)+lnr(x,y)，对取对数图像进行傅里叶变换，则

即：

Z(u,v)＝F_i(u,v)+F_r(u,v)

其中u,v为二维空频率，Z为原始图像的傅里叶变换，Fi为照度图像的傅里叶变换，Fr为反射图像的傅里叶变换。

如果借助一个滤波函数H来处理图像，那么可以较好的去除图像中的乘性噪声。

步骤2，对去噪处理后图像进行sorbel边缘提取，并使用图像形态学方法对边缘提取图像进行膨胀腐蚀操作，从而将图像填充为信号区域为凸性边缘的二值图像。

步骤3，计算整体图像的粗糙度，粗糙度代表了图像中平均纹理的大小，可以作为MeanShift算法中距离聚类阈值hs的理想选择。计算图像像素平均偏移，这个值代表图像各像素点之间像素值变化情况，可以作为MeanShift算法中像素值聚类阈值hr的理想选择。

粗糙度计算的具体方法如下。首先对每个像素点，定义其邻居窗口的平均值为M_k(x,y)，也即：

其中f(i,j)表示图像在点(i,j)的像素值，k＝1,2,……,L，其中2^k×2^k为图像最大的窗口，则L的值可由下式求出：2^L≤B＜2^L+1，B为图像的最短边的长度。

然后计算窗口不重叠的最大差值，E_k,h(x,y)，E_k,v(x,y)，E_k,d(x,y)，其中E_k,h(x,y)代表水平方向的最大差值，E_k,v(x,y)代表垂直方向的最大差值，E_k,d(x,y)代表对角方向的最大差值。这三者的计算方式如下：

E_k,h(x,y)＝|M_k(x+2^k-1,y)-M_k(x-2^k-1,y)|

E_k,v(x,y)＝|M_k(x,y+2^k-1)-M_k(x,y-2^k-1)|

E_k,d(x,y)＝|M_k(x+2^k-1,y+2^k-1)-M_k(x-2^k-1,y-2^k-1)|

然后求出对应于每个k的E的最大值，也即：

对每个像素点，求出图像区域最佳纹理大小W_best(x,y)＝2^k。最后，计算所有W_best的和的平均值，即为粗糙度F：

像素平均偏移幅度的具体计算方式如下，图像大小为M×N(像素点数)，P为某个位置像素值大小，A为整个图像的像素平均值，整个图像的像素平均偏移量

则：

步骤4，根据步骤4中计算得到的距离聚类参数hs和像素值聚类参数hr对填充图像进行MeanShift算法进行图像分块。

Mean Shift算法基本思想为假设一个d维空间，在其中任意选取一个点，以此点作为球心，作半径为h的高维球。落在球内的所有点与圆心之间都会产生向量，向量以圆心为起点。将球内所有向量相加，相加结果即为Mean shift向量。以所得的Mean Shift向量终点作为新的起点，继续之前步骤，最终向量会收敛于概率密度最大处。

对于概率密度函数f(x)，在已知n个采样点x_i，f(x)的核函数估计为：

其中h为高维球的半径，d为空间维度，w(x_i)≥0是采样点x_i的权重，K(x)是核函数，核函数需满足∫k(x)dx＝1。

定义核函数K(x)的剖面函数k(x)，即K(x)＝k(||x||²)，定义g(x)为k(x)的负导数，即g(x)＝-k'(x)，则相应的核函数G(x)＝g(||x||²)。

求概率密度函数f(x)的梯度估计即

的估计：

由g(x)＝-k'(x)，G(x)＝g(||x||²)，上式可写为：

上式第二个中括号内即为Mean Shift向量，第一个中括号内为利用核函数G(x)对f(x)的估计，因此可简写为：

其中：

这表明Mean Shift矢量的方向指向概率密度梯度的最大下降方向，即收敛点处有

对于一副图像，像素点均匀分布于图像上，不存在点的稠密性。如何定义像素点的概率密度成为关键。

定义某像素点x的概率密度方法：以x为圆心，以h为半径，对落在维度为d高维球内的点x_i制定如下规则：

(1)像素点x的值(即电磁成像中场强的模值)与像素点x_i的值越接近，定义x点概率密度越高。

(2)x_i与x位置越接近，定义x点概率密度越高。

则对于一幅图像中某点的概率密度估计式为：

其中x^L代表像素点值的大小，x^C代表像素点位置坐标，h_L为像素值高维球半径，h_C为距离高维球半径。

可以得到对某一副图像进行Mean-Shift分割的迭代公式为：

对于抛物反射面所成干扰源图像，随机抽取其中的任意一点按照上式进行MeanShift迭代，设定收敛阈值，收敛后记录其收敛点，再随机抽取剩余点，直至遍历所有点。将收敛于位置相近像素值相近的收敛点的所有点归于同一分割区域。

步骤5，根据步骤5中的分块情况，将步骤2去噪后图像不同的区域抽出信号区域，然后分块图像补充为原有图像大小，方便后面进行图像恢复。

步骤6，对分块得到的图像进行判断，如果是同一个信号分块区域则进行合并。

判断每个分块区域所在的位置，如果存在两块区域的重叠范围超过70％，那么则判断两个区域为相同的信号分割区域，对两个区域进行合并，取两个区域的最小的区域下界和最大的区域上界，得到合并的区域。

下面对发明进行仿真，以三个不同频率辐射源同时存在于一幅电磁图像中为例，其中一个源的频率为1GHz，另外两个源的频率为3GHz。距离电磁成像系统5米处的物平面上的这三个源，经过衍射受限的电磁成像系统，在像平面得到经过点扩展函数作用退化后的图像。系统采用1米×1米的成像面，成像面上分布201×201个采样点。仿真得到的电磁图像如图2所示，为无噪声情况下不同频率多辐射源电磁图像。由于仿真图像基本没有噪声，因此省略同态滤波去噪的步骤。对图像进行sobel边缘检测后图像如图3，图3为无噪声情况下不同频率多辐射源电磁图像进行sobel边缘检测图像，进行形态学处理后结果如图4，图4为图像形态学方法处理后的二值图像，进行Mean Shift聚类分区后得到三个区域，如图5，图5为经过Meanshift分割算法后得到的单信号图像，其中a为信号1提取后图像，b为信号2提取图像，c为信号3提取图像，可以看到方法比较好的对多个不同频率辐射源电磁图像进行了分区处理，为之后对各个分区进行图像重建打下基础。

接下来进行具体实验验证，以两个喇叭天线作为电磁干扰源，频率分别为1.8GHz和3GHz，使用电磁成像系统进行成像，得到的电磁图像为图6，图6电磁成像系统实验得到含噪情况下不同频率多信号源电磁图像，可以看到噪声比较大，使用同态滤波去噪然后去除图像中的直流台阶，得到的去噪图像如图7，图7为经过同态滤波后的低噪图像，然后经过形态学处理并进行Mean Shift算法分区处理，得到的分区图像如图8，图8经过Meanshift分割算法后得到的单信号图像，其中a为信号1提取后图像，b为信号2提取图像，可以看到图像分割效果比较好。然后选取相应的点扩展函数对分区区域进行LR迭代图像恢复，再将恢复后区域合并得到最终恢复区域如图9，图9为单信号图像经过LR超分辨处理，合并多个信号超分辨后图像得到最终高分辨率电磁图像，可以看到清晰的两个辐射点源。

总之，本发明当电磁探测系统得到电磁辐射图像中包含多个频率以及多个距离较远的信号区域时，如果直接对图像使用相同的点扩展函数进行图像恢复和图像重建，得到的重建图像效果较差，图像中会有比较多的噪声且一些信号会被作为噪声抑制，因此需要将图像中各个信号区域分开，对不同区域进行不同的图像恢复。对于这种分区要求，本发明首先对图像进行同态滤波抑噪，然后使用形态学方法对图像进行处理，从而使图像更易于进行Mean Shift聚类分割，然后计算图像粗糙度和像素值平均偏移作为聚类算法的阈值参数。接下了进行MS分割，并对分割后图像判断是否为同一信号区域，对于同一信号区域进行合并。对发明进行了仿真和实验验证，结果显示方法可以高效准确的对多区域多频率电磁图像进行图像抑噪并分割。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (1)

1.基于形态学方法及Meanshift算法的电磁图像分块方法，其特征在于，实现步骤如下：

步骤1，获得未知辐射源的降质图像，对图像进行对数运算，将图像转换为照度和反射两个部分相加的形式，得到取对数后图像；

步骤2，对取对数后图像进行滤波处理，滤除图像的高频噪声，从而压低亮度范围，再对图像取指数，恢复电磁功率图像，得到去噪处理后图像；

步骤3，对去噪处理后图像进行sorbel边缘提取，并使用图像形态学方法对边缘提取图像进行膨胀腐蚀操作，从而将图像填充为信号区域为凸性边缘的整体二值图像；

步骤4，计算整体二值图像的粗糙度，粗糙度代表了图像中平均纹理的大小，作为MeanShift算法中距离聚类阈值hs，然后计算图像像素平均偏移值，此偏移值代表图像各像素点之间像素值变化情况，作为MeanShift算法中像素值聚类阈值h_r；

粗糙度计算的具体方法如下，首先对每个像素点，定义其邻居窗口的平均值为M_k(x，y)，也即：

其中f(i，j)表示图像在点(i，j)的像素值，k＝1，2，......，L，其中2^k×2^k为图像最大的窗口，则L的值可由下式求出：2^L≤B＜2^L+1，B为图像的最短边的长度；

然后计算窗口不重叠的最大差值，E_k，h(x，y)，E_k，v(x，y)，E_k，d(x，y)，其中E_k，h(x，y)代表水平方向的最大差值，E_k，v(x，y)代表垂直方向的最大差值，E_k，d(x，y)代表对角方向的最大差值，这三者的计算方式如下：

E_k，h(x，y)＝|M_k(x+2^k-1，y)-M_k(x-2^k-1，y)|

E_k，v(x，y)＝|M_k(x，y+2^k-1)-M_k(x，y-2^k-1)|

E_k，d(x，y)＝|M_k(x+2^k-1，y+2^k-1)-M_k(x-2^k-1，y-2^k-1)|

然后求出对应于每个k的E的最大值，也即：

E_k＝E_max＝max(E_1，h，E_2，h......E_L，h，E_1，v，E_2，v......E_L，v，E_1，d，E_2，d......E_L，d)

对每个像素点，求出图像区域最佳纹理大小W_best(x，y)＝2^k；最后，计算所有W_best的和的平均值，即为粗糙度F：

像素平均偏移幅度的具体计算方式如下，图像大小为M×N(像素点数)，P为某个位置像素值大小，A为整个图像的像素平均值，整个图像的像素平均偏移量

则：

步骤5，根据步骤4中计算得到的距离聚类参数h_s和像素值聚类参数h_r对填充图像进行MeanShift算法进行图像分块；

步骤6，根据步骤5中的分块情况，将步骤2去噪后图像的不同的区域抽出信号区域，然后分块图像补充为原有图像大小，方便后面进行图像恢复；

步骤7，对分块得到的图像进行判断，如果是同一个信号分块区域则进行合并。

Images