CN109063564B - 一种目标变化检测方法 - Google Patents

Abstract

本发明提供一种目标变化检测方法，首先通过Gabor纹理特征提取图像显著性区域，然后通过引导滤波融合得到保留显著地物的融合影像，通过meanshift对融合影像分割得到分块，再利用HOG纹理特征计算分块纹理方差，比较得到最终变化检测结果。本发明技术方案基于图像引导滤波融合与纹理特征分析的高分辨率遥感卫星影像用于军事目标变化检测技术，它能够通过影像融合在保证变化检测精度的同时提高变化检测的效率，真正实现了军事目标变化的自动、快速、准确检测。

Description

一种目标变化检测方法

技术领域

本发明涉及一种目标变化检测方法，尤其是一种基于图像引导滤波融合与纹理特征分析的高分辨率遥感卫星影像用于军事目标变化检测的技术，属于遥感图像处理技术领域。

背景技术

作为一门正在兴起、并有着广泛应用前景的学科，遥感技术是目前为止能够提供大面积范围的动态观测数据的唯一手段，具有空间上的连续性和时间上的序列性，可以提供进行地物感知和监测的多时相图像数据。因此，分析同一地理位置的多时相遥感图像成为获取自然环境和军事目标变化信息的一种有效途径。

现已有一些利用遥感影像进行变化检测的方法被提出，包括基于像素的变化检测方法和面向对象的变化检测方法。在针对高分辨率遥感影像时，地物目标往往由一定区域大小像素构成，以像素为单位的变化检测算法没有考虑每个像素的空间邻域特征，存在一定的局限性，因此，多采用面向对象的变化检测方法。

目前面向对象的变化检测算法，总的可以分为两类，一类是对一张影像进行分割得到对象，在另一影像上搜寻该对象，这类方法只需进行一次分割，但是对分割结果精度要求很高，并且当搜索影像上有新地物产生时，并不会被检测出来；另一类是对两张影像同时进行分割，对分割结果进行比较，这类方法需要进行两次分割，并且在分割后还要进行复杂的叠置分析、区域合并等操作，且易受配准误差影响。对于军事目标的变化检测，既要求极高的精度保证信息的准确性，又要求较高的效率保证时效性，现有的两类面向对象变化检测算法并不能同时满足这两项要求。

发明内容

本发明的目的在于克服现有技术的不足，提供一种目标变化检测方法，基于图像引导滤波融合与纹理特征分析技术，利用高分辨率遥感卫星影像快速、准确、自动化地进行军事目标变化检测，它能够克服上述现有变化检测方法用于军事目标变化检测时的不足，满足高分辨率遥感卫星影像军事目标变化检测应用的需求。

本发明目的通过如下技术方案予以实现：

提供一种目标变化检测方法，包括以下步骤：

(1)获取覆盖待检测目标的两幅卫星影像A和B，两幅卫星影像A和B分别进行几何配准和辐射校正处理；

(2)分别获取卫星影像A和B的二值化显著性图pa和pb；

(3)利用二值化显著性图pa和pb对卫星影像A和B进行引导滤波，再进行加权融合，获得融合后的影像Pab；

(4)融合影像Pab进行均值漂移分割，得到m个分割图斑块S；

(5)对卫星影像A和B进行梯度方向直方图纹理特征提取，获取卫星影像A、B的纹理强度图I_A和I_B；

(6)通过纹理强度图判断m个分割图斑块S是否发生变化，输出发生变化的分割图斑块的比例。

优选的，若比例为0则待检测目标没有发生变化，否则表明待检测目标发生了变化。

优选的，卫星影像A和B的分辨率高于3m。

优选的，若卫星影像A和B分辨率不一致，通过对低分辨的卫星影像进行重采样，使得卫星影像A和B具有同样的分辨率。

优选的，步骤(2)中获取卫星影像A和B的二值化显著性图pa和pb的具体方法为：首先使用Gabor滤波算子实部对卫星影像A和B进行滤波，分别得到卫星影像A和B的Gabor强度图；再对两幅Gabor强度图进行Otsu二值化处理，得到影像A、B的二值化显著性图pa和pb。

优选的，Gabor滤波算子形式如下：

X＝x cosθ+y sinθ

Y＝-x cosθ+y sinθ

其中，(x,y)为卫星影像A或B的像素位置，(X,Y)为像素在滤波方向等效位置；λ为Gabor函数的波长；θ为Gabor函数并行条纹的方向；ψ为相位偏移；σ为高斯函数标准差；γ为空间纵横比。

优选的，二值化处理的方式为：利用t作为二值化的阈值得到二值化显著性图pa和pb，阈值t采用Otsu法确定。

优选的，卫星影像A和B进行引导滤波，采用如下公式：

q_i＝a_kI_i+b_k

其中，i是卫星影像A或B像素索引，k是窗口索引，I是引导图像的值，p_i是二值化显著性图pa或pb的第i个像素的灰度值，q_i是输出图像qa或qb的第i个像素的灰度值，a_k和b_k是当窗口中心位于k时该线性函数的系数；μ_k是I在窗口ω_k中的平均值，

是I在窗口ω_k中的方差，|ω|是窗口ω_k中像素的数量，

是二值化显著性图pa或pb在窗口ω_k中的均值，ε为滤波器滤波效果调节参数，窗口ω_k表示第k个窗口ω。

优选的，进行加权融合的方法具体为：

其中，g(i)为融合后的影像Pab在第i个像素处灰度值，g₁(i)、g₂(i)为原始影像A、B第i个像素处灰度值，w₁(i)、w₂(i)分别为影像qa，qb第i个像素处灰度值。

优选的，融合影像Pab进行均值漂移分割，得到m个分割图斑块S的具体方法为：

4.1融合影像Pab进行图像平滑滤波，对于每一个点X_i，计算联合特征空间(x,y,g)中颜色带宽和空间带宽范围内所有点的重心特征坐标X′，g为灰度值，(x,y)为Pab的像素坐标；

4.2计算重心特征坐标X′的灰度值和坐标赋给当前处理点，重复4.1进行迭代，直到达到最大迭代次数或者两次迭代获得的X′的差值小于设定阈值；

4.3进行相似区域合并；

4.4对于小于设定面积下限的区域与其相邻且最相似的区域合并。

优选的，对卫星影像A和B进行梯度方向直方图纹理特征提取，获取卫星图像A、B的纹理强度图I_A和I_B的具体方法为：

5.1利用sobel算子计算卫星图像A、B每一点x方向、y方向梯度dx、dy；

5.2计算卫星图像A、B每一点的梯度幅值M(x,y)和梯度方向θ(x,y)；

5.3计算图像A、B每一点纹理强度I(x,y):

其中，θ(x,y)、M(x,y)分别为窗口内点的梯度方向和大小，bin[τ]为梯度直方图区间，Hist[τ]为直方图每个区间内的梯度和，C为调节参数，r为方向指向。

优选的，判断m个分割图斑块S是否发生变化的具体方法为：计算卫星影像A上第m分块的纹理方差

和卫星影像B上第i分块的纹理方差

计算变化标志bChange_m：

其中，thre_SD为方差阈值。

本发明与现有技术相比具有如下优点：

(1)本发明的目标变化检测方法，具有较强的场景适应性，抗干扰能力强，运算量较小，可有效应用于遥感图像分析处理等领域。

(2)本发明利用Gabor算子对结构特征敏感的特点，通过Gabor算子提取影像上显著性地物，使得军事目标提取准确；

(3)本发明通过引导滤波使显著图上噪声得以消除，同时边缘处得以保留，消除了融合图像上的噪声，得到较为平滑的融合影像；

(4)本发明通过影像融合最大程度保留两幅影像上军事目标，在保证分割精度的基础上只需要进行一次分割，节省处理时间；

(5)本发明通过梯度直方图特征计算影像纹理复杂度，使用像素特征描述区域特征，可以更好的抵御配准误差以及光照变化带来的影响。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的二维Gabor滤波器的示意图。

图3是本发明实施例的引导滤波融合的流程图。

图4是本发明实施例用到的Meanshift分割示意图，(x,y)为像素坐标，z为像素灰度值，其中图(a)表示原始影像，图(b)表示搜寻重心的过程，图(c)表示平滑的过程，图(d)表示经过相似区域合并的结果。

图5是本发明实施例的HOG纹理特征提取的示例，其中，图(a)表示原始影像，图(b)为HOG纹理强度图。

具体实施方式

以下结合实施例和附图详细说明本发明技术方案。

本发明所采用基于图像引导滤波融合与纹理特征分析的军事目标变化检测，它能够通过影像融合在保证变化检测精度的同时提高变化检测的效率。首先采用2维Gabor滤波器对影像进行显著性区域提取。利用原始影像对显著性图进行引导滤波，并将结果作为权重图对原始影像进行加权融合。通过对融合影像meanshift分割，得到图斑。然后利用梯度直方图特征计算两张原始影像纹理强度，并比较分块纹理强度图方差，得到最终变化检测结果。

参见图1，具体实现方法包含以下步骤：

步骤1：数据准备，包括获取待进行军事目标变化检测的卫星影像，两幅分别为影像A和影像B，卫星影像分辨率应高于3m，影像覆盖区域需包含待检测军事目标，两副影像分别进行几何配准和辐射校正处理；

变化检测的数据应为同一地区不同时间获取的遥感影像，若影像分辨率不一致，可以通过重采样将其采样至同一分辨率。变化检测容易受到几何配准误差和辐射校正误差的影响，在变化检测前进行几何配准和辐射校正是必要的。

步骤2：影像显著性区域提取，包括分别使用Gabor滤波算子对影像A、B进行滤波得到Gabor强度图，然后对Gabor强度图Otsu二值化，得到影像A、B的二值化显著性图pa和pb。

二维Gabor滤波器能够很好地描述对应于空间频率(尺度)、空间位置及方向选择性的局部结构信息，常被用于纹理表示和描述。军事目标较其周围背景地物，局部结构信息往往更为丰富，通过Gabor纹理特征，可以将军事目标作为显著性目标与背景地物区分开。

二维Gabor滤波器核函数实部的数学形式为：

X＝x cosθ+y sinθ

Y＝-x cosθ+y sinθ

其中，(x,y)为原始图像的像素位置，(X,Y)为像素在滤波方向等效位置；λ为波长，以像素为单位指定，通常大于等于2，但不能大于输入图像尺寸的五分之一；θ指定了Gabor函数并行条纹的方向，它的取值为0到360度；ψ为相位偏移，取值范围为-180度到180度；σ为高斯函数标准差，通常设为0.56λ；γ为空间纵横比，决定了Gabor函数形状。具体实施时，λ取7，选取0°、30°、60°、90°、120°、150°这6个方向Gabor算子，θ为1，ψ为0，即空间纵横比1:1，相位偏移为0。

二维Gabor滤波器核函数实部的空间可视化图像如图2所示。影像A、B分别获得各自的显著图pa和pb。

计算Otsu阈值t采用如下公式，

T＝MAX{U_t0×(V_to-V_tf)²+U_t1×(V_t1-V_tf)²}

其中，U_t0表示当阈值为t时，显著像元占总像元数比例，V_t0表示显著像元Gabor强度平均值，U_t1表示非显著像元占总像元数比例，V_t1表示非显著像元平均Gabor强度平均值，V_tf表示影像平均Gabor强度，V_tf＝U_t0×V_t0+U_t1×V_t1。从最小Gabor强度值到最大Gabor强度值遍历t，当t使得类间方差T最大时，此时t作为Otsu阈值T_J；p_z表示灰度直方图，z表示灰度值。t为显著像元和非显著像元的区分阈值。

步骤3：对影像A、B的二值化显著性图pa和pb影像进行引导滤波获得引导滤波后的影像qa，qb，再进行加权融合，获得融合后的影像Pab。以原始影像作为引导图对显著性图像进行引导滤波。若直接使用显著性图作为权重图对原始影像进行加权融合，将引入大量噪声，对后续分割带来困难，降低变化检测精度。使用引导滤波的方式可以减少噪声的引入。引导滤波融合具体流程图如图3所示。

对显著性图p进行引导滤波处理时，按下式得到引导滤波后的影像q

q_i＝a_kI_i+b_k

其中，i是像素索引，k是窗口索引，I是引导图像的值，p_i是输入图像P第i个像素的灰度值，q_i是输出图像q第i个像素的灰度值，a_k和b_k是当窗口中心位于k时该线性函数的系数。μ_k是I在窗口ω_k中的平均值，

是I在窗口ω_k中的方差，|ω|是窗口ω_k中像素的数量，

是显著性图p在窗口ω_k中的均值。ε用于防止a求得的过大，也是调节滤波器滤波效果的重要参数。窗口ω_k表示第k个窗口ω。

一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，要具体求某一点的输出值时，需将所有包含该点的线性函数值平均。利用原始影像每一点灰度值、该点的线性系数，即可求得滤波后每一点灰度值。

具体实施时，根据高分辨率遥感影像的特点，引导滤波窗口ω选取7*7大小，ε值为0.01，可以得到较好的滤波效果。

按下式使用引导滤波后的显著性图作为权重图对原始影像进行加权融合，

其中，g(i)为融合后的影像在第i个像素处灰度值，g₁(i)、g₂(i)为原始影像A、B第i个像素处灰度值，w₁(i)、w₂(i)分别为影像qa，qb第i个像素处灰度值。该式的意义在于，当某像素在两景影像上都不是显著性区域时，融合影像上该像素灰度为A、B影像灰度平均值；否则根据显著性对该像素进行加权融合，特别的，当某一像素在某一期影像上显著性为零，另一影像上显著性非零，融合结果灰度值等于显著性非零像素灰度值，显著性地物得以保留。

步骤4，融合影像Pab进行均值漂移分割，得到m个分割图斑块S，包括以下子步骤：

步骤4.1，进行图像平滑滤波，滤波前的图像如图4(a)，平滑滤波后的图像如图4(b)，对于每一个点X_i，在联合特征空间(x,y,g)中计算该点颜色带宽和空间带宽范围内所有点的重心特征坐标X′，完成平滑滤波，其中(x,y)为坐标，g为灰度值，计算公式如下：

其中，X′为计算得到的更新后重心联合特征空间坐标，包括重心灰度值和坐标，X_i为窗口内每一点特征坐标为每一点的权重，X_L、X_c分别为当前重心灰度值和坐标，X_iL、X_ic为窗口内每一点的灰度值和坐标，h_L、h_C分别为颜色带宽和空间带宽，g_L、g_C为灰度和空间距离的核函数g(x)，核函数g(x)都选用均匀分布，即：

步骤4.2，如图4(c)，将计算得到重心的灰度值和坐标赋给当前处理点，采用X′的值作为当前X_L、X_c，带入公式，进行迭代，获得下一个X′。若达到最大迭代次数(默认5次)，或者两次迭代获得的X′的差值小于设定阈值，结束搜寻；否则，则转至步骤4.1；

步骤4.3，相似区域合并，在经过平滑后，如果两个区域重心之间的距离小于设定值，即则对这两个区域进行合并，这里区域的灰度值和空间坐标使用的是区域重心灰度值和空间坐标。

步骤4.4，小区域合并，根据设置的最小区域面积，在相似区域合并后，若仍有区域面积小于给定最小面积，将其合并至与其相邻且最相似区域。

步骤4.3、4.4对应图4(d)的过程。

具体实施时，颜色带宽设置为8，空间带宽设置为20，ε选取2个像素，给定最小面积阈值应小于待检测军事目标在影像上像素大小。最小面积阈值与目标大小和影像分辨率相关，目标较大，分辨率较高时，最小面积阈值可以设置较大，对于高分影像一般使用500左右。使用默认参数进行分割后，查看分割结果，若存在过分割现象，则可以调大颜色带宽、空间带宽以及最小面积阈值；若存在分割不足，则应该调小颜色带宽、空间带宽以及最小面积阈值。

步骤5，对图像A、B进行梯度方向直方图(HOG)纹理特征提取，获取图像A、B的纹理强度图I_A和I_B，如图5所示，包括以下子步骤：

首先，利用sobel算子，计算图像A、B每一点x方向、y方向梯度dx、dy，sobel算子矩阵形式如下：

x方向梯度算子为：

y方向梯度算子为：

然后，根据x、y方向的梯度算子，判断每一点的梯度方向，每30度一个方向，一共分为12个方向，具体计算式如下：

计算每一点纹理强度，具体计算公式如下；

其中，θ(x,y)、M(x,y)分别为窗口内点的梯度方向和大小，bin[τ]为梯度直方图区间，Hist[τ]为直方图每个区间内的梯度和，I(x,y)为计算得到的纹理强度，C为调节参数。具体实施时，窗口大小选取7*7，调节参数C为15。

步骤6，判断分割图斑块S是否发生变化，输出发生变化的分割图斑块的比例。对融合影像分割得到的每一分块，分别计算该块对应区域在纹理强度图上的灰度方差，通过比较每个区域的方差变化，来判定该区域是否发生了变化。具体计算式为：

其中，

为影像A上第m分块的纹理方差，

为影像B上第m分块的纹理方差，bChange_m判定第m分块是否发生变化，thre_SD为方差阈值。具体实施时，thre_SD通常设置为2.8，根据实际情况可适当调整，若影像整体辐射变化较大或者背景地物变化也较为明显，可以适当调大该阈值。如果大于阈值表明该分割图斑块发生了变化，统计发生变化的分割图斑块的比例，并将该比例输出。若比例为0则目标没有发生变化，否则表明目标发生了变化。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (10)

1.一种目标变化检测方法，其特征在于，包括以下步骤：

(1)获取覆盖待检测目标的两幅卫星影像A和B，两幅卫星影像A和B分别进行几何配准和辐射校正处理；

(2)分别获取卫星影像A和B的二值化显著性图pa和pb；

(3)利用二值化显著性图pa和pb对卫星影像A和B进行引导滤波，再进行加权融合，获得融合后的影像Pab；

(4)融合影像Pab进行均值漂移分割，得到m个分割图斑块S；

(5)对卫星影像A和B进行梯度方向直方图纹理特征提取，获取卫星影像A、B的纹理强度图I_A和I_B；

(6)通过纹理强度图判断m个分割图斑块S是否发生变化，输出发生变化的分割图斑块的比例；

判断m个分割图斑块S是否发生变化的具体方法为：计算卫星影像A上第m分块的纹理方差

和卫星影像B上第i分块的纹理方差

计算变化标志bChange_m：

其中，thre_SD为方差阈值；如果大于方差阈值表明该分割图斑块发生了变化，统计发生变化的分割图斑块的比例，并将该比例输出；若比例为0则目标没有发生变化，否则表明目标发生了变化。

2.如权利要求1所述的目标变化检测方法，其特征在于，卫星影像A和B的分辨率高于3m。

3.如权利要求2所述的目标变化检测方法，其特征在于，若卫星影像A和B分辨率不一致，通过对低分辨的卫星影像进行重采样，使得卫星影像A和B具有同样的分辨率。

4.如权利要求1所述的目标变化检测方法，其特征在于，步骤(2)中获取卫星影像A和B的二值化显著性图pa和pb的具体方法为：首先使用Gabor滤波算子实部对卫星影像A和B进行滤波，分别得到卫星影像A和B的Gabor强度图；再对两幅Gabor强度图进行Otsu二值化处理，得到卫星影像A和B的二值化显著性图pa和pb。

5.如权利要求4所述的目标变化检测方法，其特征在于，Gabor滤波算子形式如下：

X＝x cosθ+y sinθ

Y＝-x cosθ+y sinθ

其中，(x，y)为卫星影像A或B的像素位置，(X，Y)为像素在滤波方向等效位置；λ为Gabor函数的波长；θ为Gabor函数并行条纹的方向；ψ为相位偏移；σ为高斯函数标准差；γ为空间纵横比。

6.如权利要求4或5所述的目标变化检测方法，其特征在于，二值化处理的方式为：利用t作为二值化的阈值得到二值化显著性图pa和pb，阈值t采用Otsu法确定。

7.如权利要求1或2所述的目标变化检测方法，其特征在于，卫星影像A和B进行引导滤波，采用如下公式：

q_i＝a_kI_i+b_k

其中，i是卫星影像A或B像素索引，k是窗口索引，I是引导图像的值，p_i是二值化显著性图pa或pb的第i个像素的灰度值，q_i是输出图像qa或qb的第i个像素的灰度值，a_k和b_k是当窗口中心位于k时该线性函数的系数；μ_k是I在窗口ω_k中的平均值，

是I在窗口ω_k中的方差，|ω|是窗口ω_k中像素的数量，

是二值化显著性图pa或pb在窗口ω_k中的均值，ε为滤波器滤波效果调节参数，窗口ω_k表示第k个窗口ω。

8.如权利要求7所述的目标变化检测方法，其特征在于，进行加权融合的方法具体为：

其中，g(i)为融合后的影像Pab在第i个像素处灰度值，g₁(i)、g₂(i)为原始影像A、B第i个像素处灰度值，w₁(i)、w₂(i)分别为影像qa，qb第i个像素处灰度值。

9.如权利要求8所述的目标变化检测方法，其特征在于，融合影像Pab进行均值漂移分割，得到m个分割图斑块S的具体方法为：

4.1融合影像Pab进行图像平滑滤波，对于每一个点X_i，计算联合特征空间(x，y，g)中颜色带宽和空间带宽范围内所有点的重心特征坐标X′，g为灰度值，(x，y)为Pab的像素坐标；

4.2将重心特征坐标X′的灰度值和坐标赋给当前处理点，重复4.1进行迭代，直到达到最大迭代次数或者两次迭代获得的X′的差值小于设定阈值；

4.3进行相似区域合并；

4.4对于小于设定面积下限的区域与其相邻且最相似的区域合并。

10.如权利要求9所述的目标变化检测方法，其特征在于，对卫星影像A和B进行梯度方向直方图纹理特征提取，获取卫星图像A、B的纹理强度图I_A和I_B的具体方法为：

5.1利用sobel算子计算卫星图像A、B每一点x方向、y方向梯度dx、dy；

5.2计算卫星图像A、B每一点的梯度幅值M(x，y)和梯度方向θ(x，y)；

5.3计算图像A、B每一点纹理强度I(x，y)：

其中，θ(x，y)、M(x，y)分别为窗口内点的梯度方向和大小，bin[τ]为梯度直方图区间，Hist[τ]为直方图每个区间内的梯度和，C为调节参数，τ为方向索引。

Images