本发明属于遥感图像处理及图像识别技术领域，具体涉及一种基于联合显著性分析的sar图像油库区检测方法。
背景技术
合成孔径雷达(syntheticapertureradar，sar)不受天气、光照等条件的限制，具有无可比拟的全天候、全天时等诸多优点，因此被广泛应用于船只检测、地形勘察、资源勘探、目标识别以及环境监测等各个方面。sar系统中雷达技术的不断发展使得sar图像的分辨率日益提高，数据量不断增加，因此传统的sar图像解译方法已经不能满足海量数据的处理要求。如何从复杂的地物场景中有效检测出sar图像目标，如：车辆、舰船、坦克、油库等，目前已成为遥感图像处理及图像识别技术领域的研究热点。
最早开展sar图像目标识别技术研究的是美国麻省理工学院的林肯实验室。林肯实验室的saip系统，使用的是经典的模板匹配方法。基于模板匹配的方法简单易实现，但由于sar图像特征常常受各种因素的影响发生显著的变化，使得模板匹配方法在面对较为复杂的目标或有较强背景杂波干扰时常常发生漏检与误检等问题，已无法满足人们的需要。
一些专家学者借鉴了光学遥感影像的目标检测方法，提出了sar图像的目标检测模型。但是sar图像成像机制与光学遥感影像不同，如果直接使用现有的模型来处理sar图像，将存在以下几个主要问题：(1)由于sar独特的相干成像原理，sar图像中存在大量的斑点噪声，这和光学遥感影像中的加性噪声不同。斑点噪声严重地影响了sar图像的质量，并且也给后续的目标检测带来不利影响。(2)sar数据成像时对目标的方位角变化很敏感，方位角不同会导致成像结果产生很大的差异，sar图像中的目标信息也会产生变化。因此光学遥感影像中常用的边界信息无法在sar图像中发挥作用。(3)与光学遥感影像不同，sar图像缺少光谱信息，因此一些常用的光学遥感影像特征无法应用于sar图像。此外，sar图像中的目标一般在图像中只占很小的比例，使得目标检测精度严重下降。例如，使用zhang等人在文章“regionsofinterestdetectioninpanchromaticremotesensingimagesbasedonmultiscalefeaturefusion”中提出的基于多尺度特征融合的mff模型进行sar图像目标检测，将无法获取精确结果。
人眼在面对外界海量的信息时处理能力是有限的，视觉注意机制是保证信息处理过程中高效率的关键，它使得人眼只对其中显著的信息做出响应而舍弃其他不重要的信息。视觉显著性就是受此原理发展起来的，它的目标是在处理图像时使计算机可以具备类似人眼视觉的主动性和选择性。由于显著性目标检测具有不依赖先验知识、识别准确率高以及能快速检测出图像中的显著目标等优势，它已经被广泛应用于自然景物图像目标检测和光学遥感影像分析领域。利用显著性检测方法获取图像中的显著目标后，其结果可用于图像重建、图像融合和变化检测等领域，并且可以大大降低计算量。例如imamoglu等人在文章“asaliencydetectionmodelusinglow-levelfeaturesbasedonwavelettransform”中提出的基于小波变换的wt模型，能够较为准确高效的计算自然景物图像的显著图并提取显著目标。因此，将视觉显著性引入sar图像处理，能够为高效、精确地从sar图像中提取出目标提供新的思路。
但是，将视觉显著性引入sar图像目标检测还有一些问题需要注意。一方面sar图像中缺少由可见光提供的真实色彩信息，因此需要针对sar图像寻找新的特征，另一方面一组具有相似特征的sar图像间存在共性关系，结合多幅图像中的共性显著特征进行目标检测，不仅保留了视觉显著性分析的优点，而且能够有效降低由真实色彩缺失所带来的sar图像目标检测精度降低问题，还能够进一步抑制背景干扰，避免误判。联合显著性分析就是指在一组具有相似特征的图像中挖掘共有的显著目标，因此本发明使用联合显著性分析对sar图像进行目标检测以提高检测精度。
技术实现要素：
本发明的目的在于提供一种基于联合显著性分析的sar图像油库区检测方法，该方法用于对一组sar图像中的油库区进行精确检测。传统的sar图像目标检测主要借鉴光学遥感影像的目标检测方法，首先需要对sar图像中的斑点噪声进行建模，但是对sar图像斑点噪声的建模是一个很复杂的问题，需要考虑多方面因素构建模型。此外sar图像的成像机理不同于光学遥感影像，针对sar图像寻找合适的特征也是一个重要问题。因此本发明方法主要关注三个方面：
1)降低斑点噪声对检测结果的不利影响；
2)针对sar图像缺少光谱信息问题，寻找合适的新的特征；
3)挖掘多幅图像间的共有信息弥补sar图像颜色信息的缺失，实现更为准确的sar图像的油库区检测。
本方法针对一组具有相似地物特征的sar图像进行处理，首先，利用方向平滑滤波器对一组sar图像进行降噪处理，其次，提取一组sar图像的亮度特征、纹理特征和曲线特征，进行共性显著特征分析，得到一组sar图像的共性显著特征图，再次，基于共生直方图，计算得到一组sar图像中每一幅图像的单图显著特征图，然后，融合共性显著特征图和单图显著特征图，得到一组sar图像的最终显著图，最后通过最大类间方差法对最终显著图进行阈值分割提取油库区，具体包括以下步骤：
步骤一：进行sar图像降噪处理，即构建一个增强的方向平滑滤波器，通过该滤波器的核函数在滑动窗口中执行低通滤波，使一组sar图像中每幅图像的斑点噪声均得到有效降低；
步骤二：在一组sar图像中，提取每幅图像各像素的亮度特征、纹理特征和曲线特征，即对于一组sar图像中的每一幅图像，提取图像中各像素的亮度信息作为该像素的亮度特征，对图像进行下采样，在下采样图像上使用基于邻域阈值的局部二值模式提取每个像素的纹理特征，对图像中的所有像素点进行张量投票，计算得到每个像素的曲线特征，将提取的各像素的亮度特征、纹理特征以及曲线特征作为本发明中聚类使用的共性显著特征；
步骤三：利用模糊c均值聚类算法完成聚类，即在一组sar图像中，利用每幅图像各像素的亮度特征、纹理特征和曲线特征，构建像素的共性显著特征向量，通过模糊c均值聚类算法，将一组sar图像中所有像素的共性显著特征向量进行聚类，得到k个簇；
步骤四：计算一组sar图像中每幅图像的共性显著特征图，即通过测量各像素的共性显著特征向量在各个簇中的出现频率，来计算簇与簇之间的距离，然后将第i个簇中含有的像素数与总像素数相除，相除的结果定义为第i个簇的权重，其中，得到所有k个簇的权重后，令一个簇到其他簇的距离的加权和为该簇的对比度显著值，把簇的对比度显著值赋给每一个属于该簇的像素点，由此获得一组共性显著特征图；
步骤五：计算一组sar图像中每幅图像的单图显著特征图，即在一组sar图像的每一幅图像中，统计所有强度值为a的像素的8邻域内，强度值为b的像素的个数，其中a和b的取值范围为0至255，然后将统计得到的像素对的个数排列成一个方阵，该方阵即为图像的共生直方图，然后用共生直方图中不同数值的出现概率的负对数，来描述直方图中每个点的初始显著值，最后增强初始显著值差异，得到一组sar图像的单图显著特征图；
步骤六：计算一组sar图像中每幅图像的最终显著图，即对一组sar图像的每一幅单图显著特征图进行归一化，将归一化后的单图显著特征图中每个像素的灰度值记为对应sar图像相同位置像素的显著性权重，用该权重来增强对应sar图像的共性显著特征图，从而获得一组sar图像的最终显著图；
步骤七：利用最大类间方差法准确提取一组sar图像中每幅图像的油库区，即通过最大类间方差法计算得到最终显著图的分割阈值，利用该阈值将一组sar图像中每幅图像的最终显著图分割为一幅二值图像模板，该模板中“1”代表油库区，“0”代表非油库区，最后将二值图像模板与对应sar图像相乘，从而得到该sar图像的油库区提取结果。
附图说明
图1为本发明的流程图。
图2为本发明的sar图像共性显著特征聚类结果。第一行是sar图像，第二行是聚类结果。
图3为本发明的共性显著性分析与单图显著性分析的结果对比。第一行是共性显著特征图sp_m；第二行是单图显著特征图sp_s。
图4为示例sar图像及其对应的地面实况(ground-truth)。第一行是示例sar图像；第二行是地面实况图。
图5为示例图片采用本发明方法和其他方法生成的显著图的比较。(a)-(d)分别表示采用mff，wt，csfa和本发明方法的结果，其中csfa是发明内容中步骤二至步骤四所提出的共性显著特征分析方法(commonsalientfeatureanalysis,csfa)。
图6为示例图片采用本发明方法和其他方法所检测到的油库区比较。(a)-(d)分别表示采用mff，wt，csfa和本发明方法的结果。
具体实施方式
下面结合附图对本发明做进一步详细说明。本发明的总体框架如图1所示，现介绍每一步的实现细节。
步骤一：进行sar图像降噪处理；
本发明中使用的一组尺寸为m×n的sar图像，sar图像中的斑点噪声通常被建模为纯乘性噪声过程，用下式表示
v(r,c)＝l(r,c)s(r,c)
其中v表示雷达测量到的值，l表示图像的实际辐射值，s表示斑点噪声。对于单视sar图像，s是瑞利分布或均值为1的负指数分布。对于具有独立外观的多视sar图像，s是均值为1的伽马分布。m表示图像的行数，n表示图像的列数，r、c分别表示像素的横、纵坐标，r＝1,2，...,m，c＝1,2，...,n。
增强的方向平滑(enhanceddirectionalsmoothing,eds)滤波器可以降低图像中的斑点噪声，它用核函数在滑动窗口中执行低通滤波。滤波窗口的典型尺寸范围为3×3到33×33，窗口的大小通常设为奇数。较大的滤波器窗口意味着可以利用较大的图像区域进行计算，但可能需要更多的计算时间，具体取决于滤波器算法的复杂程度。如果滤波器窗口的尺寸过大，重要细节将因过度平滑而丢失。另一方面，如果滤波器窗口的尺寸太小，则不利于降低斑点噪声。根据相关文献的实验结果分析，实验中采用3×3或7×7的滤波窗通常会产生最佳结果。本发明中使用的是一组尺寸为m×n的sar图像，利用尺寸为3×3的滤波器降噪，将降噪后的sar图像组记为
步骤二：提取一组sar图像中每幅图像所有像素的亮度特征、纹理特征和曲线特征；
在本方法中，将滤波后的sar图像的亮度特征、纹理特征以及油罐特有的曲线特征作为本发明中聚类使用的共性显著特征。
亮度特征：直接使用ix中像素pn的亮度信息作为该像素的亮度特征gn。pn表示图像ix中的第n个像素，n＝1，2，…,m×n。
纹理特征：纹理特征是从ix中提取的，主要包括两步。在本发明中，首先对ix以1/8的比例进行下采样，然后在下采样后的图像上使用基于邻域阈值的局部二值模式来获取纹理特征图，计算每个像素的纹理特征texturen的方法如下：
texturen＝∑ss(|gnw-gn|-δ)
其中gn和gnw分别表示当前像素及其邻域像素的灰度值。本发明中，pnw表示当前像素8-领域内的像素，δ是一个阈值，将其设为40。最后将纹理特征图恢复至初始大小就可获得每个像素的纹理特征。
曲线特征：张量投票是一种可以从噪声图像中提取显著特征结构的算法，因此本发明利用张量投票来获取sar图像中每个像素的曲线特征curven。二维平面中张量投票的投票域有两种：棒形张量和球形张量。
首先，对图像ix中的所有像素点进行球形张量编码，即对其编码为：
根据矩阵谱分解原理，张量tensor可分解如下的特征值λ和对应的特征向量e的线性组合。
tensor＝(λ1-λ2)e1e1t+λ2(e1e1t+e2e2t)
其中，λ1和λ2为张量的非负特征值；e1和e2为两个特征值对应的特征向量；e1e1t和e1e1t+e2e2t分別表示棒形张量和球形张量；对应的(λ1-λ2)和λ2分别表示线特征的显著性和点特征的显著性。因此，利用下式即可获得每个像素的曲线特征curven。
curven＝λ1-λ2
步骤三：利用模糊c均值聚类算法完成聚类；
经过步骤二后，我们得到了gn，texturen和curven三种sar图像的共性显著特征。考虑到每个特征对聚类结果的贡献值，本发明对每个特征赋予不同的权重，然后将其组合为一个显著特征向量featuren。
featuren＝(0.8·gn,0.8·texturen,curven)
模糊c均值聚类(fuzzyc-meansalgorithm，fcm)是一种柔性划分方法，主要是度量每个数据点属于某个簇的隶属度。fcm这种对样本类别划分的不确定性是一种更为客观的描述，因此本发明采用fcm完成多幅sar图像的共性显著特征聚类。fcm使用u函数来确定每个数据属于某个簇的程度。
其中k表示簇的个数。ui(pn)表示像素pn属于第i个簇ci的隶属度。fcm的目标函数j如下：
公式中的m是加权指数，默认m＝2。ci是簇ci的中心。fcm的算法步骤如下：首先，初始化各个簇中心ci使其满足
步骤四：计算一组sar图像中每幅图像的共性显著特征图；
全局对比度的计算基于一个通用规则：具有独特特征的稀有群集在直观上是引人注目的并且是更加显著的。为了达到这个目的，本发明首先构造了一个特征空间，空间的三个维度分别对应三个共性显著特征：gn，texturen和curven。对每个维度分别进行均匀量化，量化后的级数分别表示为bg，bt和bc，其中，本发明设置bg＝bt＝bc＝8。
对于每个簇，通过测量各特征向量在该簇中的出现频率来计算它在特征空间中的特征直方图。基于特征直方图，簇ci的对比度显著值s(ci)可以被估计为在特征空间中到其他簇的距离的加权和，具体计算如下。
ω(ci)表示簇ci中的像素数量与所有图像的总像素数量之比。d(ci,cj)表示簇ci与簇cj的距离，ci表示第i个簇，bn代表直方图中特征量化的级数，bn＝bg×bt×bc。另外，fi,t表示第t个特征向量出现在簇ci中的频率。由此可获得共性显著特征图sp_m。
步骤五：计算一组sar图像中每幅图像的单图显著特征图；
共生直方图是表示像素之间的共生关系的一种统计方法；它可以揭示图像灰度值的空间关系的本质。共生直方图中的元素是图像中像素的强度信息。如果图像是由灰度值相似的像素块组成，则共生直方图的主对角线上元素将是相对较大的值；反之如果图像像素灰度值在局部有波动，则共生直方图的非主对角线元素具有相对较大的值。
对于一幅m×n的图像，共生直方图coh计算如下：
coh＝[coh(a,b)]
coh是一个256×256的对称矩阵。coh(a,b)就表示强度值为a的像素的8-邻域内，强度值为b的像素的个数，强度值的取值范围为0至255。
coh就是由所有的coh(a,b)构建而成，然后使用下面的公式将coh归一化到0-1。
pcoh(a,b)表示的是共生直方图中不同数值的出现概率。理论上，对于研究者想要的显著地感兴趣区域，这些区域应该具有一个较小的pcoh(a,b)，相反，对于非显著的区域，pcoh(a,b)应该是较大的数值。
为了使pcoh(a,b)与强度对(a,b)的显著性的关系是一个正相关的映射，受boltzmann熵定理的启发，用负对数关系来描述每个亮度对的显著性关系，如下公式所示：
l(a,b)＝-ln(pcoh(a,b))
其中，l(a,b)可以被看作是图像的初始显著值。因此显著区域应该具有一个较大的l(a,b)，与此同时，对于非显著区域而言l(a,b)应该是一个较小的数值。但是，此时显著区和非显著区的初始显著值的差别不是很大，因此考虑增强这两者之间的差异。本发明利用k-means算法来增强两者之间的差异：首先，利用k-means将两者区分开，然后分别对两类做相应计算来增加两者之间的差异。增强计算如下：
其中，es(a,b)表示的是最后增强后的显著值。y(a,b)代表k-means算法的分类结果，值为1说明在该像素在较大数值的簇里，值为0则说明该像素位于较小数值的类别中。u表示共生直方图的平均显著值，经过此操作，显著区域与非显著区域之间的差异就会被进一步加大。
单图显著特征图中每个像素的显著值都取决于它的邻域，对于像素pn，它的显著值smn计算如下：
smn＝∑es(gn,gnw)
由此得到的显著图中包含一些道路和其他冗余信息。因此，本发明用高斯滤波器来处理显著图，在一定程度上削弱了这些非显著区域的强度，由此可获得单图显著特征图sp_s。
步骤六：计算一组sar图像中每幅图像的最终显著图；
从图3可以看出：第一，单图显著性分析方法在显著背景干扰区域上产生了严重的误判，而共性显著性分析则准确地抑制其显著性。第二，相比于sp_m，sp_s中背景干扰更加明显。第三，虽然sp_m能有效地抑制背景干扰，但是显著性较高的油库区中的部分油罐并不完整。
综上，这两种方法的检测结果都有各自的优点，所以本发明考虑融合这两种显著特征图来去除所有多余的信息的干扰。为了避免油库区中的部分油罐不完整，本发明使用单图显著特征图sp_s来优化共性显著特征图sp_m，以计算最终显著图msp。
msp＝γ×sp_m
γ＝sp_s/255
第一个公式的作用是融合单图显著特征图和共性显著特征图获得最终显著图，这里是一种点乘运算。首先对sp_s进行归一化获得单幅图像中内部显著性的权重指标γ，然后用γ来增强sp_m获得最终显著图msp，如图5所示。从图中可以看出，msp确实明显地突显出油库区，有效地抑制了非油库区的显著性。
步骤七：利用最大类间方差法准确提取sar图像的油库区；
通过最大类间方差法得到最终显著图的分割阈值，利用该阈值将最终显著图分割为一幅二值图像模板，用“1”代表油库区，用“0”代表非油库区。最后将二值图像模板与sar图像进行点乘得到最终的油库区提取结果。
本发明的效果可通过以下实验结果与分析进一步说明：
1.实验数据
实验中用的sar数据是2014年8月29日在东京地区由alos2卫星拍摄获取的。该sar图像是由卫星的l波段采用hh极化方式获得，分辨率为3m。我们从实验数据中随机选取7幅图像进行结果展示，图4展示了示例sar图像及其对应的地面实况(ground-truth)。
2.对比实验
为了评价本发明方法的性能，我们设计了如下的对比实验，选取了现有的具有代表性的视觉注意方法选取了多特征融合方法(mff)，小波变换方法(wt)，以及步骤二至步骤四所提出的共性显著性分析方法(csfa)与本发明方法进行性能对比。从主观上分别对比了不同方法生成的显著图和感兴趣区域图，如图5和图6所示。
从结果可以看出：除了csfa和本发明方法，其他的方法在显著背景干扰区域上都产生了严重的误判。mff，wt的显著性检测结果中都包含有大量的背景信息，而csfa和本发明则对背景信息进行了有效地抑制。这是由于后者采用对所有显著特征进行聚类然后计算显著值，可以有效地将油库区和背景区域进行区分。而前四种方法主要针对单图中的显著特征进行计算，常用的亮度以及频域分析无法准确地区分与油库区表现相似的背景区域。