CN101950364A - 基于邻域相似度和阈值分割的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于邻域相似度和阈值分割的遥感图像目标变化检测方法，主要解决传统方法对含有较高噪声的目标变化检测抗噪性能差，检测精度低的缺点。其实现过程包括：(1)利用强度归一化公式对两幅遥感图像进行灰度匹配；(2)利用邻域相似性距离测度构造两幅遥感图像的相似矩阵；(3)结合相似矩阵构造两幅遥感图像的差异影像；(4)对差异影像构造二维灰度直方图，采用2D-OTSU方法确定分割阈值，将目标区域和背景区域分离；(5)用模糊熵的方法对未处理的边缘和噪声点继续进行分类。本发明具有抗噪性强，检测变化目标精确度高的优点，可用于检测多时相遥感图像变化的目标。

Description

基于邻域相似度和阈值分割的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及图像变化检测方法，适用于两幅同一地区不同时相遥感图像中含有地物变化信息的检测。

背景技术

变化检测旨在通过同一地区不同时期的图像间的差异来得到图像的地物变化信息，它是开展森林资源调查、土地利用、覆盖变化研究、环境灾害评估、城市规划和国防军情监控等对地观测应用中的关键技术，具有迫切的需要和广泛的应用前景。

变化检测是遥感领域研究的重点和热点问题，许多学者从不同的角度对现有的变化检测方法进行了分类、分析。变化检测方法主要考虑图像对图像的检测，是在像元级水平上发展起来的检测方法。现有的变化检测方法可以归纳为两大类：一类是监督检测法，另一类是非监督检测法。前者是指根据地面真实数据来获取变化区域的训练样区，从而进行变化检测；后者是直接对两个不同时相的数据检测而不需要任何额外的信息。由于地面的真实信息不容易得到，因此非监督的变化检测方法是常用的变化检测方法。非监督检测法通常做法是直接比较同一位置不同时相的像元特征值来检测变化，通常采用数学变换的方式产生不同时相间的差异影像，再对差异影像进行阈值化处理，从中提取变化区域。在现阶段，从数据源的角度还提出了影像和数字线划图、影像和地图、DEM、DOM与影像的变化检测方法。从影像是否配准角度出发，提出了先配准影像再变化检测和变化检测与影像配准同时进行的方法。从算法的角度而言，许多学者提出了多种方法的综合使用，这些方法包括使用人工神经网络、马尔科夫随机场、数学形态学等。

尽管各种变化检测方法已经广泛的应用到诸多领域，如EM-MRF、OTSU、聚类等方法。然而，从整体上来说变化检测仍然存在着不少的困难和问题有待深入的研究和解决。如差异影像的构建，在较强噪声的情况下，传统的滤波方法容易造成模糊效应，丢失目标信息；在对图像的分类处理中，上述方法的抗噪性不够强，无法在目标边缘保留完好的同时有效的去除噪声干扰，或者容易造成目标边界信息的丢失，导致检测结果的效果不佳，检测精度不够高。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于邻域相似度和阈值分割的遥感图像变化检测方法，以在保留变化区域的同时有效的去除噪声干扰，提高目标的检测精度。

实现本发明目的技术方案是先利用像素的邻域相似性度量公式构造差异图，然后使用2D-OTSU方法进行分割，最后用模糊熵对分割的结果进行后处理，得到最后结果，其具体步骤包括如下：

(1)选取不同时相的两幅遥感图像X₁和X₂，时相中包含有目标变化信息，且目标信息受到噪声污染；

(2)采用强度归一化法对两幅遥感图像进行灰度匹配；

(3)对灰度匹配后的两幅遥感图像分别取同一位置像素点的局部邻域块，设定窗口尺度，用滑窗法对每个像素点取相同尺度邻域，对于每个邻域，度量两幅图像对应位置像素点的相似性，得到两幅图像的相似矩阵：

d_a(i，j)＝∫_ΩG_a(s)|X(x+s)-Y(y+s)|²ds

其中X和Y分别为两幅遥感图像相同位置对应的邻域块，s为窗口尺度，x和y为窗口中心像素点的位置，i，j为相似矩阵中像素点坐标，Ω表示窗口区域，G_a(s)为高斯核函数；

(4)利用相似矩阵定义每个像素点的权值：

$w(i,j)=1-e^{{-d}_{\mathrm{\α}}(i,j)/h}$

其中h为衰减参数；

利用定义的权值，求取两幅遥感图像的差异影像：

X_D(i，j)＝w(i，j)×(X₁(i，j)-X₂(i，j))

X₁(i，j)和X₂(i，j)分别为两幅遥感图像在位置(i，j)的像素灰度；

(5)设定任意阈值为(s，t)，用2D-Otsu方法求得差异影像的最佳阈值(s₀，t₀)，根据最佳阈值将背景和目标分割开来，即将s≥s₀，t≥t₀的点作为目标，将s＜s₀，t＜t₀的点作为背景，将s＞s₀，t≤t₀以及s≤s₀，t＞t₀的点作为边缘点和噪声点；

(6)使用模糊熵处理方法对边缘点和噪声点继续分割，检测出目标。

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

1.本发明由于利用像素点的邻域块的相似性度量构造相似矩阵，利用相似矩阵定义的权值并同时结合图像差值，使构造得到的差异影像能够在保留目标区域比较完好的同时，极大的减小了噪声的影响，相比传统的滤波处理效果更好；

2.本发明由于采用2D-OTSU方法对图像进行分割，考虑到了图像的空间信息和像素点邻域之间的相关性，因而相比传统的单阈值分割方法，在信噪比较低的情况下，能够减小噪声对分割结果的影响，提高目标的检测精度；

3.本发明由于在阈值分割后，采用模糊熵的方法对未分割的边缘点和噪声点进行分类，弥补忽略边界信息所带来的误差，更具有合理性。

附图说明

图1是本发明的实现流程示意图；

图2是本发明实验使用的两幅不同时相遥感图像；

图3是本发明实验选取的图像变化区域参考图像；

图4是本发明实验用距离度量公式构造的相似矩阵图；

图5是本发明实验利用相似矩阵构造的差异图；

图6是本发明实验的差异影像的二维灰度直方图；

图7是用2D-OTSU方法对差异影像分割得到的初始分割图；

图8是用模糊熵对初始分割图进行再分类得到的最终结果图。

具体实施方式

参照图1，本发明的具体实施步骤如下：

步骤1.选取两幅不同时相的遥感图像；

实验图像选取的是意大利Elba岛西部地区的两幅不同时相的遥感图像，如图2所示，其中图2(a)为第一时相的遥感图像，图2(b)为第二时相的遥感图像，这两幅不同时相的遥感图像中包含有目标变化信息，且目标信息受到噪声污染。

步骤2.用强度归一化公式对两幅遥感图像进行灰度匹配；

采用强度归一化法对两幅遥感图像进行灰度匹配，即通过下式对两幅遥感图像进行强度归一化：

X₂＝(σ₁/σ₂)×(X₂-μ₂)+μ₁

其中μ₁，μ₂分别为两幅遥感图像X₁，X₂的均值，σ₁，σ₂分别为两幅遥感图像X₁，X₂的方差。

步骤3.构造两幅遥感影像的相似矩阵

设定窗口的尺度，分别取两幅遥感图像同一位置的局部邻域块，用滑窗法对每个像素点取相同尺度邻域，对于每个邻域，度量两幅图像对应位置像素点的相似性，得到两幅图像的相似矩阵：

d_α(i，j)＝∫_ΩG_α(s)|X(x+s)-Y(y+s)|²ds

其中X和Y分别为两幅遥感图像相同位置对应的邻域块，x和y为窗口中心像素点的位置，i，j为相似矩阵中像素点坐标，Ω表示窗口区域，G_α(s)为高斯核函数，s为窗口尺度，取值为5×5，相似矩阵用图像形式表示如图4，从图4中可以看出，相似矩阵能够有效的突出变化目标区域，很大程度削弱了噪声的干扰。

步骤4.利用相似矩阵构造差异影像。

利用相似矩阵定义权值：

$w(i,j)=1-e^{{-d}_{\mathrm{\α}}(i,j)/h}$

其中h为衰减参数，实验中h＝400，利用权值以及两图像的差值，求取两幅遥感图像的差异影像：

X_D(i，j)＝w(i，j)×(X₁(i，j)-X₂(i，j))

差异影像如图5所示，从图5所示的差异影像效果，更能看出，噪声得到了很大的去除，同时变化目标也得到了很好的保留。

步骤5.用2D-OTSU方法对差异影像进行分割。

5a)构造差异影像的二维灰度直方图

设差异影像的大小为m×n，灰度级L＝256，差异图像中任一像素点的坐标为(i，j)，其灰度值为f(i，j)，1≤i≤m，1≤j≤n，定义像素点(i，j)的邻域平均灰度为g(i，j)，g(i，j)的计算公式如下：

$g(i,j)=\frac{1}{D\×D}\underset{x^{\′}=-(D-1)/2}{\overset{(D-1)/2}{\mathrm{\Σ}}}\underset{y^{\′}=-(D-1)/2}{\overset{(D-1)/2}{\mathrm{\Σ}}}f(i+x^{\′},j+y^{\′})$

其中，D为邻域窗口宽度，实验中窗口D设为5×5大小，

由差异图像灰度值f(i，j)和其邻域灰度值g(i，j)同时出现的次数，可构成二维灰度直方图L_i，j，0≤i，j≤L-1，如图6所示。

5b)用2D-OTSU方法对差异影像进行分割

5b1)将二维灰度直方图L_i，j归一化为概率p_ij＝L_i，j/(m×n)；

5b2)设定(s，t)为任意门限向量，0≤s，t＜L，利用二维灰度直方图的归一化概率p_ij，分别计算背景区域和目标区域发生的概率分布为：

${\mathrm{\ω}}_0(s,t)=\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{p}_{\mathrm{ij}},$

${\mathrm{\ω}}_1(s,t)=\underset{i=s+1}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}},$

其中ω₀(s，i)表示背景区域概率分布，ω₁(s，t)表示目标区域概率分布；

5b3)计算背景均值矢量μ₀和目标对应的均值矢量μ₁分别为：

${\mathrm{\μ}}_0={[{\mathrm{\μ}}_{0i},{\mathrm{\μ}}_{0j}]}^T={[\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}/{\mathrm{\ω}}_0(s,t),\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}/{\mathrm{\ω}}_0(s,t)]}^T$

${\mathrm{\μ}}_1={[{\mathrm{\μ}}_{1i},{\mathrm{\μ}}_{1j}]}^T={[\underset{i=s+1}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}/{\mathrm{\ω}}_1(s,t),\underset{i=s+1}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}/{\mathrm{\ω}}_1(s,t)]}^T$

5b4)根据归一化为概率p_ij计算总体均值矢量μ_z为：

${\mathrm{\μ}}_z={[{\mathrm{\μ}}_{\mathrm{zi}},{\mathrm{\μ}}_{\mathrm{zj}}]}^T={[\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}},\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}]}^T$

5b5)定义一个目标和背景类间的离散测度矩阵：

σ_B＝ω₀(s，t)[(μ₀-μ_z)(μ₀-μ_z)^T]+ω₁(s，t)[(μ₁-μ_z)(μ₁-μ_z)^T]则采用矩阵σ_B的迹trσ_B(s，t)作为目标和背景类间的距离测度函数：

$\mathrm{tr}{\mathrm{\σ}}_B(s,t)=\frac{[({\mathrm{\ω}}_0(s,t){\mathrm{\μ}}_{\mathrm{zi}}-{\mathrm{\μ}}_i(s,t){)}^2+({\mathrm{\ω}}_0(s,t){\mathrm{\μ}}_{\mathrm{zj}}-{\mathrm{\μ}}_j(s,t){)}^2]}{{\mathrm{\ω}}_0(s,t)(1-{\mathrm{\ω}}_0(s,t))}$

其中， ${\mathrm{\μ}}_{\mathrm{zi}}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}},$ ${\mathrm{\μ}}_{\mathrm{zj}}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}},$ ${\mathrm{\μ}}_i(s,t)=\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}},$ ${\mathrm{\μ}}_j(s,t)=\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}j{p}_{\mathrm{ij}},$

5b6)将在trσ_B(s，t)取最大值时的阈值，作为2D-OTSU图像分割法的最佳阈值(s₀，t₀)，即：

trσ_B(s₀，t₀)＝MAX{trσ_B(s，t)}，0≤s，t≤L

5b7)将s≤s₀，t≤t₀的点定为背景点，将s≥s₀，t≥t₀的点定为目标点，将s＞s₀，t≤t₀以及s≤s₀，t＞t₀的点定为边缘点或者噪声点，得到的初始分割结果如图7所示。

步骤6.用模糊熵的方法对初始分割后的图像进行后处理

对初始分割后的图像进行后处理是对除背景点和目标点以外的边缘点和噪声点进行再分割，具体步骤如下：

6a)设Q＝[q(x，y)]_M×N为2D-OTSU初步分割过的图像，取一个大小为n×n，n为奇数，中心在(x，y)的窗口W_n(x，y)，实验中窗口大小取为3×3，令m₀＝0代表背景，m₁＝255代表目标，用如下公式分别表示窗口中像素对背景区域和目标区域的隶属程度：

${\mathrm{\μ}}_{m_0}\left(q(x+k,y+l)\right)={[1+|q(x+k,y+l)-m_0|/\mathrm{\λ}]}^{-1}$

${\mathrm{\μ}}_{m_1}\left(q(x+k,y+l)\right)={[1+|q(x+k,y+l)-m_1|/\mathrm{\λ}]}^{-1}$

其中

q(x+k，y+l)为在(x+k，y+l)位置处像素的灰度值；λ为一常数，在本实验中取值为255；

6b)定义背景模糊熵和目标模糊熵分别为：

$e\left(m_0\right)=\frac{1}{(n\×n)\ln 2}\underset{k=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}\underset{l=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}S\left\{{\mathrm{\μ}}_{m_0}\right[q(x+k,y+l)\left]\right\}$

$e\left(m_1\right)=\frac{1}{(n\×n)\ln 2}\underset{k=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}\underset{l=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}S\left\{{\mathrm{\μ}}_{m_1}\right[q(x+k,y+l)\left]\right\}$

其中 $S\left\{{\mathrm{\μ}}_{m_i}\right[q(x+k,y+l)\left]\right\}={-\mathrm{\μ}}_{m_i}\left[q(x+k,y+l)\right]\×\ln \left\{{\mathrm{\μ}}_{m_i}\right[q(x+k,y+l)\left]\right\}$

$-\{1-{\mathrm{\μ}}_{m_i}\[q(x+k,y+l)]\}\×\ln \{1-{\mathrm{\μ}}_{m_i}\left[q(x+k,y+l)\right]\},i=\mathrm{0,1}$

对同一个窗口W_n(x，y)中的像素，求出e(m₀)和e(m₁)并比较两者大小，如果最小值为e(m₀)，则令q(x，y)为背景m₀，如果最小值为e(m₁)，则令q(x，y)为目标m₁，得到最终的分割图像如图8所示。

本发明的效果可以通过以下实验进一步说明。

1)实验条件

本实验采用Landsat 5ETM拍摄的意大利Elba岛西部地区1994年8月和1994年9月的两幅大小为326×414的遥感图像作为实验数据，采用matlab 7.0软件作为仿真工具。计算机配置为Intel Core2/1.8G/1G。

2)实验内容和检测结果

本实验先通过构造实验数据中的两幅遥感图像的差异影像，再对差异影像分割成目标区域和背景区域，得到的目标区域为变化类，背景区域为非变化类，变化检测实验结果如图8所示，从图8可见，本实验的方法在大量去除噪声的同时，使目标的边缘信息得到较好的保留，检测结果在视觉效果上良好。

3)检测精度评价

选用如图3所示的二值图作为参考图，将本实验的检测结果与参考图进行对比，得出检测精度评价如表1。

表1检测精度评价

	错误个数	正检率	虚警率
				Elba岛西部地区	295	92.71％	5.05％

从表1可见，本实验方法的错误个数较少，正检率高，虚警率低，检测精度较高。

Claims (4)

1.一种基于邻域相似度和阈值分割的遥感图像变化检测方法，包括如下步骤：

(1)选取不同时相的两幅遥感图像X₁和X₂，时相中包含有目标变化信息，且目标信息受到噪声污染；

(2)采用强度归一化法对两幅遥感图像进行灰度匹配；

(3)对灰度匹配后的两幅遥感图像分别取同一位置像素点的局部邻域块，设定窗口尺度，用滑窗法对每个像素点取相同尺度邻域，对于每个邻域，度量两幅图像对应位置像素点的相似性，得到两幅图像的相似矩阵：

d_a(i，j)＝∫_ΩG_a(s)|X(x+s)-Y(y+s)|²ds

其中X和Y分别为两幅遥感图像相同位置对应的邻域块，s为窗口尺度，x和y为窗口中心像素点的位置，i，j为相似矩阵中像素点坐标，Ω表示窗口区域，G_a(s)为高斯核函数；

(4)利用相似矩阵定义每个像素点的权值：

其中h为衰减参数；

利用定义的权值，求取两幅遥感图像的差异影像：

X_D(i，j)＝w(i，j)×(X₁(i，j)-X₂(i，j))

其中X₁(i，j)和X₂(i，j)分别为两幅遥感图像在位置(i，j)的像素灰度；

(5)设定任意阈值为(s，t)，用2D-Otsu方法求得差异影像的最佳阈值(s₀，t₀)，根据最佳阈值将背景和目标分割开来，即将s≥s₀，t≥t₀的点作为目标，将s＜s₀，t＜t₀的点作为背景，将s＞s₀，t≤t₀以及s≤s₀，t＞t₀的点作为边缘点和噪声点；

(6)使用模糊熵处理方法对边缘点和噪声点继续分割，检测出目标。

2.根据权利要求1所述的图像目标变化检测方法，其中步骤(2)所述的采用强度归一化法对两幅遥感图像进行灰度匹配，是通过下式进行强度归一化：

X₂＝(σ₁/σ₂)×(X₂-μ₂)+μ₁

其中μ₁，μ₂分别为两幅遥感图像X₁，X₂的均值，σ₁，σ₂分别为两幅遥感图像X₁，X₂的方差。

3.根据权利要求1所述的图像目标变化检测方法，其中步骤(5)所述的用2D-Otsu方法求得最佳阈值(s₀，t₀)，按如下步骤进行：

3a)设(s，t)为任意门限向量，0≤s，t＜L-1，灰度级数L＝256，求解图像进行分割时目标区域和背景区域发生的概率分别为：

${\mathrm{\ω}}_0(s,t)=\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{p}_{\mathrm{ij}},$ ${\mathrm{\ω}}_1(s,t)=\underset{i=s+1}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}$

其中p_ij为各灰度级的归一化概率；

3b)计算背景均值矢量μ₀和目标对应的均值矢量μ₁分别为：

${\mathrm{\μ}}_0={[{\mathrm{\μ}}_{0i},{\mathrm{\μ}}_{0j}]}^T={[\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}/{\mathrm{\ω}}_0(s,t),\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}/{\mathrm{\ω}}_0(s,t)]}^T$

${\mathrm{\μ}}_1={[{\mathrm{\μ}}_{1i},{\mathrm{\μ}}_{1j}]}^T={[\underset{i=s+1}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}/{\mathrm{\ω}}_1(s,t),\underset{i=s+1}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}/{\mathrm{\ω}}_1(s,t)]}^T$

3c)根据归一化概率p_ij计算总体均值μ_z为：

${\mathrm{\μ}}_z={[{\mathrm{\μ}}_{\mathrm{zi}},{\mathrm{\μ}}_{\mathrm{zj}}]}^T={[\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}},\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}]}^T$

3d)定义一个目标和背景类间的离散测度矩阵：

σ_B＝ω₀(s，t)[(μ₀-μ_z)(μ₀-μ_z)^T]+ω₁(s，t)[(μ₁-μ_z)(μ₁-μ_z)^T]，

采用离散测度矩阵σ_B的迹trσ_B(s，t)作为目标和背景类间的距离测度函数：

$\mathrm{tr}{\mathrm{\σ}}_B(s,t)=\frac{[({\mathrm{\ω}}_0(s,t){\mathrm{\μ}}_{\mathrm{zi}}-{\mathrm{\μ}}_i(s,t){)}^2+({\mathrm{\ω}}_0(s,t){\mathrm{\μ}}_{\mathrm{zj}}-{\mathrm{\μ}}_j(s,t){)}^2]}{{\mathrm{\ω}}_0(s,t)(1-{\mathrm{\ω}}_0(s,t))}$

其中 ${\mathrm{\μ}}_{\mathrm{zi}}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}},$ ${\mathrm{\μ}}_{\mathrm{zj}}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}j{p}_{\mathrm{ij}},$ ${\mathrm{\μ}}_i(s,t)=\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}},$ ${\mathrm{\μ}}_j(s,t)=\underset{i=0}{\overset{s}{\mathrm{\Σ}}}\underset{j=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}};$

3e)将在trσ_B(s，t)取最大值时的阈值，作为2D-OTSU图像分割法的最佳阈值(s₀，t₀)，即：

trσ_B(s₀，t₀)＝MAX{trσ_B(s，t)}，0≤s，t＜L。

4.根据权利要求1所述的图像目标变化检测方法，其中步骤(6)所述的使用模糊熵处理方法对边缘点和噪声点继续分割，按如下步骤进行：

4a)设Q＝[q(x，y)]_M×N为2D-OTSU初步分割过的图像，取一个大小为n×n，n为奇数，中心在(x，y)的窗口W_n(x，y)，令m₀＝0代表背景，m₁＝255代表目标，用如下公式分别表示窗口中像素对背景区域和目标区域的隶属程度：

${\mathrm{\μ}}_{m_0}\left(q(x+k,y+l)\right)={[1+|q(x+k,y+l)-m_0|/\mathrm{\λ}]}^{-1}$

${\mathrm{\μ}}_{m_1}\left(q(x+k,y+l)\right)={[1+|q(x+k,y+l)-m_1|/\mathrm{\λ}]}^{-1}$

其中

q(x+k，y+l)为在(x+k，y+l)位置处像素的灰度值；λ为一常数，在本实验中取值为255；

4b)定义背景模糊熵和目标模糊熵分别为：

$e\left(m_0\right)=\frac{1}{(n\×n)\ln 2}\underset{k=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}\underset{l=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}S\left\{{\mathrm{\μ}}_{m_0}\right[q(x+k,y+l)\left]\right\}$

$e\left(m_1\right)=\frac{1}{(n\×n)\ln 2}\underset{k=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}\underset{l=-\frac{n-1}{2}}{\overset{\frac{n-1}{2}}{\mathrm{\Σ}}}S\left\{{\mathrm{\μ}}_{m_1}\right[q(x+k,y+l)\left]\right\}$

其中 $S\left\{{\mathrm{\μ}}_{m_i}\right[q(x+k,y+l)\left]\right\}={-\mathrm{\μ}}_{m_i}\left[q(x+k,y+l)\right]\×\ln \left\{{\mathrm{\μ}}_{m_i}\right[q(x+k,y+l)\left]\right\}$

$-\{1-{\mathrm{\μ}}_{m_i}\[q(x+k,y+l)]\}\×\ln \{1-{\mathrm{\μ}}_{m_i}\left[q(x+k,y+l)\right]\},i=\mathrm{0,1}$

对同一个窗口W_n(x，y)中的像素，求出e(m₀)和e(m₁)并比较两者大小，如果最小值为e(m₀)，则令q(x，y)为背景m₀，如果最小值为e(m₁)，则令q(x，y)为目标m₁。

Images