CN101738607A - 基于聚类的高阶累量交叉熵的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于聚类的高阶累量交叉熵的SAR图像变化检测方法，主要解决现有模型法检测中，在变化范围很小的情况下检测结果不准确的缺点。其实现过程包括：(1)利用K均值聚类方法提取窗口的像素点；(2)利用欧氏距离去掉多余像素点；(3)用提取出的窗口的像素点作为样本点进行基于高阶累量的交叉熵算法计算，得到变化差异图；(4)对差异图用最大类间方差阈值方法得到二值图；(5)对二值图像考虑邻域关系进行基于邻域的去杂点运算得到最终的变化检测结果。本发明具有检测小目标变化的优点，可用于检测多时相SAR图像变化的目标。

Description

基于聚类的高阶累量交叉熵的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及图像目标变化检测方法，适用于两幅不同时相SAR图像中含有目标变化信息，被噪声严重污染的变化检测。

背景技术

变化检测旨在通过同一地区不同时期的图像间的差异来得到感兴趣的地物变化信息，它是开展森林资源调查、土地利用、覆盖变化研究、环境灾害评估、城市规划和国防军情监控等对地观测应用中的关键技术，具有迫切的需要和广泛的应用前景。合成孔径雷达SAR具有全天候，全天时的特点，是很好的变化检测信息源，研究SAR图像变化检测技术有着非常广阔的应用前景。

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

尽管各种变化检测方法已经广泛的应用到诸多领域，然而，从整体上来说变化检测对小目标的检测效果仍然存在着变化细节严重丢失和在变化范围很小的情况下检测结果不准确的问题和二值图像中目标无法明确从杂点中提取的问题。

发明内容

本发明的目的在于克服上述已有问题的缺点，提出了一种基于聚类的高阶累量交叉熵的SAR图像变化检测方法，以自动实现SAR图像对小目标的变化检测，并提出了一种在二值图像中将目标从杂点中提取出来的方法。

实现本发明目的技术方案是首先分别对两幅不同时相的图像取对应窗口，分别检测两幅图像中对应窗口中与窗口中心点属同一类别的像素，然后通过两幅图像对应窗口像素分布类型的差异得到图像的变化差异图，用最大类间方差阈值方法对变化差异图进行分割，得到包含目标的二值图像。最后对二值图像中每一个像素点进行判断，得到变化目标像素点。其具体步骤包括如下：

(1)选取两幅不同时相的SAR图像，对单幅图像以每一个像素为中心选取窗口作为该中心点的邻域，对窗口内像素进行K均值聚类，根据聚类结果进行类别合并，确定窗口内类别，提取与中心像素点同属一类的像素点；

(2)将步骤(1)中提取出的第一幅图的像素点个数n₁和第二幅图的像素点个数n₂进行个数一致性处理，如果n₁＜n₂，则在第二幅图窗口中计算所有与窗口中心点同类的像素点与窗口中心点的欧式距离，并选取n₁个欧式距离最小的像素点作为第二幅图窗口的像素点；如果n₁＞n₂，则在第一幅图窗口中计算所有与窗口中心点同类的像素点与窗口中心点的欧式距离，并选取n₂个欧式距离最小的像素点作为第一幅图窗口的像素点，以处理后的两组像素点作为样本点进行基于高阶累量的交叉熵计算，得到变化差异图；

(3)对每一个像素点重复步骤(1)～(2)，得到一幅变化差异图；

(4)对得到的变化差异图用最大类间方差阈值方法取阈值，得到包含目标像素点的二值图像；

(5)在二值图像中以每一个像素为中心选取窗口边长为L的正方形窗口，窗口内像素个数为L²；

(6)计算窗口内部包含的目标像素点的个数N，如果N＞L+1，则确认此窗口中心像素是目标像素点；如果N≤L+1，计算与窗口中心像素点相连通的目标像素点的个数n，如果n＞(L+1)/2，则确认此窗口中心像素点是目标像素点，如果n≤(L+1)/2，则将此窗口中心像素点标记为被错误检测的像素点，将所有被标记的像素点的像素值设为0；

(7)重复步骤(5)～(6)，对每一个像素点进行判断，得到最终的变化检测结果图。

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

1.由于基于模型分布差异的变化检测方法对像素分布类型依赖性较大，而小尺度图像与大尺度图像相比在相同尺度范围内所有像素点都隶属于同一分布类型的置信度很小，如果窗口内所有像素点都参与变化检测，则在检测过程中极易产生细节丢失的情况，本发明利用K均值聚类提取窗口内相同分布类型的像素点，用提取出的像素点参与计算，极大地提高了检测效果；

2.本发明采用K均值聚类选择像素点后，对区域块再单独作变化检测，提高了检测精度；

3.由于SAR图像含有高噪声，在小目标变化检测中容易错误的将干扰像素当成小目标像素检测出来，本发明利用变化目标像素存在连通性的特点，提出了去除干扰像素点的方法，提高了检测精度。

附图说明

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

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

图3是本发明与现有对比实验方法的差异图像；

图4是通过本发明得到的含有杂点的变化检测二值图像；

图5是通过本发明得到的遥感图像变化检测结果与现有对比实验方法的变化检测结果图以及变化检测结果参考图。

具体实施方式

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

步骤1.选取两幅不同时相的SAR图像。

实验图像选取的是两幅不同时相的水域目标SAR图像，如图2所示，其中图2(a)为第一时相，图2(b)为第二时相，变化信息为水域目标的变化，这两种不同时相的SAR图像中包含有目标变化信息，且两幅SAR图像的大小均为318×318，记图像边长为R，R＝318。

步骤2.在窗口内提取与窗口中心点同类别的像素点，具体步骤如下：

1a)以坐标(i，j)为中心，分别取边长为R/6的窗口W₁和W₂，用K均值聚类方法将窗口内像素聚为M类，M∈[6，10]；

1b)分别计算各个类别的像素值均值，并按均值大小对各个类别进行排序，将排序后的各类别的均值记为μ_i(i＝1，2，…，M)；

1c)计算两个相邻类别的均值之差Δμ_ij((i，j)＝(1，2)，(2，3)，…，(M-1，M))和均值之差的均值μ；

1d)比较相邻两个类别的均值之差Δμ_ij((i，j)＝(1，2)，(2，3)，…，(M-1，M))与均值之差的均值μ，如果第i类和第j类的均值之差Δμ_ij小于0.8μ或大于1.2μ则将i和j合并为一类，否则不合并；

1e)在两个窗口分别取与对应中心点同属一类的所有像素点集合，记为J₁和J₂，其包含的像素个数分别为n₁和n₂。

步骤3.将提取出的第一幅图像像素点集合J₁的像素点个数n₁和第二幅图像像素点集合J₂的像素点个数n₂进行个数一致性处理，处理后的像素点集合分别记为X和Y。

该一致性处理，根据第一个窗口的像素点个数n₁和第二个窗口的像素点个数n₂的大小不同，处理过程如下：

3a)如果n₁＜n₂，则计算第二幅图像窗口中的像素点集合J₂中的所有像素点与窗口中心点的欧式距离d，假设像素点I∈J₂，I的坐标为(i，j)，窗口中心像素点C的坐标为(m，n)，则像素点I和窗口中心点的欧氏距离 $d=\sqrt{{(i-m)}^2+{(j-n)}^2},$ 在集合J₂中选取n₁个欧式距离最小的像素点作为第二幅图窗口W₂的像素点，此时将两幅图像窗口中选取得到的像素点集合分别记为X和Y；

3b)如果n₁＞n₂，则计算第一幅图像窗口中的像素点集合J₁中的所有像素点与窗口中心点的欧式距离d，假设像素点J∈J₁，J的坐标为(i，j)，窗口中心像素点C的坐标为(m，n)，则像素点J和窗口中心点的欧氏距离 $d=\sqrt{{(i-m)}^2+{(j-n)}^2},$ 在集合J₁中选取n₂个欧式距离最小的像素点作为第一幅图窗口W₁的像素点，此时将两幅图像窗口中选取得到的像素点集合分别记为X和Y。

步骤4.用在两幅图像窗口中选取得到的像素点集合X和Y中的像素点进行基于高阶累量的交叉熵计算，得到变化差异图，如图3(a)所示，具体步骤如下：

4a)计算集合X的一阶矩m_X；1和二阶矩m_X；2，并根据m_X；1和m_X；2计算K_X；1和K_X；2，方程式如下：

$m_{X;1}={\∫}_0^{255}x{f}_X\left(x\right)\mathrm{dx}$

$m_{X;2}={\∫}_0^{255}{x}^2{f}_X\left(x\right)\mathrm{dx},$

K_X；1＝m_X；1

$K_{X;2}=m_{X;2}-m_{X;1}^2$

f_X(x)为集合X的概率密度函数，x∈[0，255]；

4b)计算集合Y的一阶矩m_Y；1和二阶矩m_Y；2，并根据m_Y；1和m_Y；2计算K_Y；1和K_Y；2，方程式如下：

$m_{Y;1}={\∫}_0^{255}x{f}_Y\left(x\right)\mathrm{dx}$

$m_{Y;2}={\∫}_0^{255}{x}^2{f}_Y\left(x\right)\mathrm{dx},$

K_Y；1＝m_Y；1

$K_{Y;2}=m_{Y;2}-m_{Y;1}^2$

f_Y(x)为集合Y的概率密度函数，x∈[0，255]；

4c)根据集合X、K_X；1和K_X；2计算集合X的变形X′，根据集合Y、K_Y；1和K_Y；2计算集合Y的变形Y′：

$X^{\′}=(X-K_{X;1})K_{X;2}^{-1/2}$

                   ；

$Y^{\′}=(Y-K_{Y;1})K_{Y;2}^{-1/2}$

4d)计算集合X′的n阶矩m_X′；n(n＝1，2，3，4)，并根据m_X′；n(n＝1，2，3，4)计算K_X′；3和K_X′；4，方程式如下：

$m_{X^{\′};1}={\∫}_0^{255}x{f}_{X^{\′}}\left(x\right)\mathrm{dx}$

$m_{X^{\′};2}={\∫}_0^{255}{x}^2{f}_{X^{\′}}\left(x\right)\mathrm{dx}$

$m_{X^{\′};3}={\∫}_0^{255}{x}^3{f}_{X^{\′}}\left(x\right)\mathrm{dx},$

$m_{X^{\′};4}={\∫}_0^{255}{x}^4{f}_{X^{\′}}\left(x\right)\mathrm{dx}$

$K_{X^{\′};3}=m_{X^{\′};3}-3m_{X^{\′};2}{m}_{X^{\′};1}+2m_{X^{\′};1}^3$

$K_{X^{\′};4}=m_{X^{\′};4}-4m_{X^{\′};3}{m}_{X^{\′};1}-3m_{X^{\′};2}^2+12m_{X^{\′};2}{m}_{X^{\′};1}^2-6m_{X^{\′};1}^4$

f_X′(x)为集合X′的概率密度函数，x∈[0，255]；

4e)计算集合Y′的n阶矩m_Y′；n(n＝1，2，3，4)，并根据m_Y′；n(n＝1，2，3，4)计算K_Y′；3和K_Y′；4，方程式如下：

$m_{Y^{\′};1}={\∫}_0^{255}x{f}_{Y^{\′}}\left(x\right)\mathrm{dx}$

$m_{Y^{\′};2}={\∫}_0^{255}{x}^2{f}_{Y^{\′}}\left(x\right)\mathrm{dx}$

$m_{Y^{\′};3}={\∫}_0^{255}{x}^3{f}_{Y^{\′}}\left(x\right)\mathrm{dx},$

$m_{Y^{\′};4}={\∫}_0^{255}{x}^4{f}_{Y^{\′}}\left(x\right)\mathrm{dx}$

$K_{Y^{\′};3}=m_{Y^{\′};3}-3m_{Y^{\′};2}{m}_{Y^{\′};1}+2m_{Y^{\′};1}^3$

$K_{Y^{\′};4}=m_{Y^{\′};4}-4m_{Y^{\′};3}{m}_{Y^{\′};1}-3m_{Y^{\′};2}^2+12m_{Y^{\′};2}{m}_{Y^{\′};1}^2-6m_{Y^{\′};1}^4$

f_Y′(x)为集合Y′的概率密度函数，x∈[0，255]；

4f)计算n阶切比雪夫多项式H_n(x)(n＝3，4，6)，方程式如下：

${(-1)}^n\frac{d^n\mathrm{\α}\left(x\right)}{{\mathrm{dx}}^n}=H_n\left(x\right)\·\mathrm{\α}\left(x\right)$

$\mathrm{\α}\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{x^2}{2}},$

其中，x∈[0，255]；

4g)计算含有n个像素的集合X的均值μ_X和方差σ² _X：

${\mathrm{\μ}}_X=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i}{n}$ (x_i∈X|i＝1，2，…，n)，

${\mathrm{\σ}}_X^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-{\mathrm{\μ}}_X)}^2$ (x_i∈X|i＝1，2，…，n)

计算含有m个像素的集合Y的均值和方差μ_Y和σ² _Y：

${\mathrm{\μ}}_Y=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{y}_i}{m}$ (y_i∈Y|i＝1，2，…，n)；

${\mathrm{\σ}}_Y^2=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(y_i-{\mathrm{\μ}}_Y)}^2$ (y_i∈Y|i＝1，2，…，n)

4h)用均值μ_X、方差σ² _X计算高斯函数G_X(x)，用均值μ_Y、方差σ² _Y计算高斯函数G_Y(x)：

$G_X\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_X}{e}^{-\frac{{(x-{\mathrm{\μ}}_X)}^2}{2{\mathrm{\σ}}_X^2}}$

$G_Y\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_Y}{e}^{-\frac{{(x-{\mathrm{\μ}}_Y)}^2}{2{\mathrm{\σ}}_Y^2}},$

其中，x∈[0，255]；

4i)计算集合X的概率密度f_X(x)和集合Y的概率密度f_Y(x)：

$f_X\left(x\right)=(1+\frac{K_{X^{\′};3}}{6}{H}_3\left(x\right)+\frac{K_{X^{\′};4}}{24}{H}_4\left(x\right)+\frac{10K_{X^{\′};3}^2}{720}{H}_6\left(x\right))G_X\left(x\right)$

$f_Y\left(x\right)=(1+\frac{K_{Y^{\′};3}}{6}{H}_3\left(x\right)+\frac{K_{Y^{\′};4}}{24}{H}_4\left(x\right)+\frac{10K_{Y^{\′};3}^2}{720}{H}_6\left(x\right))G_Y\left(x\right),$

其中，x∈[0，255]；

4j)用集合X的概率密度f_X(x)和集合Y的概率密度f_Y(x)计算交叉熵：

$K\left(Y\right|X)=\∫\log \frac{f_X\left(x\right)}{f_Y\left(x\right)}{f}_X\left(x\right)\mathrm{dx}$

$K\left(X\right|Y)=\∫\log \frac{f_Y\left(x\right)}{f_X\left(x\right)}{f}_Y\left(x\right)\mathrm{dx},$

KLD＝K(Y|X)+K(X|Y)

其中x∈[0，255]，

KLD为该窗口中心点在差异图像中的像素值；

4h)计算所有像素点的KLD，得到一幅变化差异图像。

步骤5.计算差值图的最大类间方差阈值，根据此阈值将差值图中感兴趣区域中大于阈值的像素点作为变化类，否则作为非变化类，得到一幅二值图像，该二值图像包含真实目标与伪目标，如图4所示。

步骤6.对二值图像考虑邻域关系进行基于邻域的去杂点运算得到最终的变化检测结果，具体步骤如下：

6a)在二值图像中以每一个像素为中心选取窗口边长为L的正方形窗口，窗口内像素个数为L²；

6b)计算窗口内部包含的目标像素点的个数N；

6c)根据N的大小采用不同的处理方法：

如果N＞L+1，则确认此窗口中心像素是目标像素点；

如果N≤L+1，由于在目标邻域有连通性的情况下，有极大地置信度可以认为该目标是真实目标，而在目标邻域基本没有连通性的情况下，可以认为该目标是伪目标，所以计算与窗口中心像素点相连通的目标像素点的个数n，如果n＞(L+1)/2，则确认此窗口中心像素点是目标像素点；如果n≤(L+1)/2，则将此窗口中心像素点标记为被错误检测的像素点，将所有被标记的像素点的像素值设为0；

6d)对每一个像素点进行判断，得到最终的变化检测结果图，如图5(a)所示，从图5(a)的结果可以看出本发明比较精确的检测出了变化目标。

本发明的效果通过以下实验具体说明：

1.实验数据

本发明使用了两组遥感图像。第一组真实实验数据集是通过Landsat 7ETM+4波段获得的墨西哥郊外的两幅大小为512×512的真实图像，256灰度级，第一时相图像是在2000年4月获得，如图2(a)所示；第二时相SAR图像是在2002年5月获得，如图2(b)所示，变化区域主要为大火破坏了大面积的当地植被所致。图5(c)为图2(a)和图2(b)变化检测结果参考图，图5(a)为本发明关于图2(a)和图2(b)的变化检测结果。第二组数据是常用于遥感变化检测精度评价的模拟MTRSI数据集。该数据集由ATM(Airborne Thematic Mapper)3波段影像和模拟变化影像构成，其中，ATM影像位于英国Feltwell村庄的一个农田区，模拟变化影像通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到，影像大小均为470×335，256灰度级，两幅影像配准误差为1.5个像元左右。

2.对比实验及实验评价指标

Jordi Inglada等学者在2007年文章“A NEW Statistical Similarity Measurefor Change Detection in Multitemporal SAR Images and Its Extension toMultiscale Change Analysis”中提出一种基于高阶累量交叉熵变化检测方法，该方法将两时相遥感图像进行运算得到差异图。

本发明设计一个实验来验证本发明的有效性。为了验证在得到变化差异图的方面本发明的有效性，实验采用本文检测变化差异的方法和基于高阶累量的交叉熵的检测变化差异的方法进行对比，来验证本发明检测出的感兴趣区域对变化检测结果的影响，其中本发明的变化差异图，如图3(a)所示，基于高阶累量的交叉熵方法得到的变化差异图，如图3(b)所示，本发明的变化检测结果图，如图5(a)所示，基于高阶累量的交叉熵方法得到的变化检测结果图，如图5(b)所示。

实验中的对变化检测结果进行量和质的分析。量的评价指标包括虚警数、漏检数以及总错误数，质的评价是将变化检测结果图5(a)、图5(b)与参考图5(c)进行主观视觉对比。

3.实验结果与分析

实验中，本发明得到的变化检测差异图如图3(a)所示，基于高阶累量交叉熵变化检测方法得到的变化检测差异图如图3(b)所示。从图3(a)与图3(b)的比较可以看出，图3(a)中可以明显看出变化目标的大概范围，而图3(b)中几乎检测不出变化目标。本发明得到的变化检测结果如图5(a)所示，基于高阶累量交叉熵变化检测方法得到的变化检测结果如图5(b)所示。从图5(a)中可以看出，本发明很接近于参考图5(c)，而图5(b)中变化信息缺失极大。

表1给出了实验中第一组实验数据结果技术指标，表2给出了实验中第二组实验数据结果技术指标。从表1和表2可以看出，对两组真实的遥感图像，本发明的总错误数比基于高阶累量交叉熵方法的总错误数分别少了66％和82％，明显的检测出变化区域。说明变化区的完整性得到了保持。总的来说，本发明实现了对小目标的变化检测，同时较好的在去除杂点的过程中保持了变化的细节。

表1实验的第一组实验数据结果

表2实验的第二组实验数据结果

Claims (5)

1.一种基于聚类的高阶累量交叉熵的SAR图像变化检测方法，包括如下步骤：

(1)选取两幅不同时相的SAR图像，对单幅图像以每一个像素为中心选取窗口作为该中心点的邻域，对窗口内像素进行K均值聚类，根据聚类结果进行类别合并，确定窗口内类别，提取与中心像素点同属一类的像素点；

(2)将步骤(1)中提取出的第一幅图的像素点个数n₁和第二幅图的像素点个数n₂进行个数一致性处理，如果n₁＜n₂，则在第二幅图窗口中计算所有与窗口中心点同类的像素点与窗口中心点的欧式距离，并选取n₁个欧式距离最小的像素点作为第二幅图窗口的像素点；如果n₁＞n₂，则在第一幅图窗口中计算所有与窗口中心点同类的像素点与窗口中心点的欧式距离，并选取n₂个欧式距离最小的像素点作为第一幅图窗口的像素点，以处理后的两组像素点作为样本点进行基于高阶累量的交叉熵计算，得到变化差异图；

(3)对每一个像素点重复步骤(1)～(2)，得到一幅变化差异图；

(4)对得到的变化差异图用最大类间方差阈值方法取阈值，得到包含目标像素点的二值图像；

(5)在二值图像中以每一个像素为中心选取窗口边长为L的正方形窗口，窗口内像素个数为L²；

(6)计算窗口内部包含的目标像素点的个数N，如果N＞L+1，则确认此窗口中心像素是目标像素点；如果N≤L+1，计算与窗口中心像素点相连通的目标像素点的个数n，如果n＞(L+1)/2，则确认此窗口中心像素点是目标像素点，如果n≤(L+1)/2，则将此窗口中心像素点标记为被错误检测的像素点，将所有被标记的像素点的像素值设为0；

(7)重复步骤(5)～(6)，对每一个像素点进行判断，得到最终的变化检测结果图。

2.根据权利要求1所述的基于聚类的高阶累量交叉熵的SAR图像变化检测方法，其中步骤(1)所述的对窗口内像素进行K均值聚类，是先取原SAR图像边长的1/6作为窗口大小；再采用K均值聚类方法，将窗口内像素聚为M个类别，M∈[6，10]。

3.根据权利要求1所述的基于聚类的高阶累量交叉熵的SAR图像变化检测方法，其中步骤(1)所述的对窗口内像素聚类的结果进行类别合并，按如下步骤进行：

3a)计算各类内像素的均值，按均值大小将各类排序，将排序后的各类别均值记为μ_i(i＝1，2，…，M)；

3b)计算相邻的两个类别的均值之差Δμ_ij((i，j)＝(1，2)，(2，3)，…，(M-1，M))和全部相邻类别的均值之差的均值μ；

3c)比较相邻的两个类别的均值之差Δμ_ij((i，j)＝(1，2)，(2，3)，…，(M-1，M))与均值之差的均值μ，如果第i类和第j类的均值之差Δμ_ij小于0.8μ或大于1.2μ，则将i和j合并为一类，否则不合并。

4.根据权利要求1所述的基于聚类的高阶累量交叉熵的SAR图像变化检测方法，其中步骤(2)所述的像素点(i，j)与像素点(m，n)之间的欧式距离d用方程式表示如下：

d²＝(i-m)²+(j-n)²    (4A)

5.根据权利要求1所述的基于聚类的高阶累量交叉熵的SAR图像变化检测方法，其中步骤(2)所述的基于高阶累量的交叉熵计算得到变化差异图，按如下步骤进行：

5a)计算两个对应窗口像素集合X的概率密度f_X和集合Y的概率密度f_Y：

$f_X\left(x\right)=(1+\frac{K_{X^{\′};3}}{6}{H}_3\left(x\right)+\frac{K_{X^{\′};4}}{24}{H}_4\left(x\right)+\frac{10K_{X^{\′};3}^2}{720}{H}_6\left(x\right))G_X\left(x\right)---\left(5A\right)$

$f_Y\left(x\right)=(1+\frac{K_{Y^{\′};3}}{6}{H}_3\left(x\right)+\frac{K_{Y^{\′};4}}{24}{H}_4\left(x\right)+\frac{10K_{Y^{\′};3}^2}{720}{H}_6\left(x\right))G_Y\left(x\right)---\left(5B\right)$

其中，K_X；n为集合X的n阶矩m_X；n组成的表达式：

K_X；1＝m_X；1

$K_{X;2}=m_{X;2}-m_{X;1}^2---\left(5C\right)$

$K_{X;3}=m_{X;3}-{3m}_{X;2}{m}_{X;1}+2m_{X;1}^3$

$K_{X;4}=m_{X;4}-{4m}_{X;3}{m}_{X;1}-3m_{X;2}^2+12m_{X;2}{m}_{X;1}^2-6m_{X;1}^4$

K_Y；n为集合Y的n阶矩m_Y；n组成的表达式：

K_Y；1＝m_Y；1

$K_{Y;2}=m_{Y;2}-m_{Y;1}^2---\left(5D\right)$

$K_{Y;3}=m_{Y;3}-{3m}_{Y;2}{m}_{Y;1}+2m_{Y;1}^3$

$K_{Y;4}=m_{Y;4}-{4m}_{Y;3}{m}_{Y;1}-3m_{Y;2}^2+12m_{Y;2}{m}_{Y;1}^2-6m_{Y;1}^4$

$X^{\′}=(X-K_{X;1})K_{X;2}^{-1/2}---\left(5E\right)$

$Y^{\′}=(Y-K_{Y;1})K_{Y;2}^{-1/2}---\left(5F\right)$

H_n(x)是n阶切比雪夫多项式，由下式计算：

${(-1)}^n\frac{d^n\mathrm{\α}\left(x\right)}{{\mathrm{dx}}^n}=H_n\left(x\right)\·\mathrm{\α}\left(x\right)---\left(5G\right)$

$\mathrm{\α}\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{x^2}{2}}$

G_X(x)和G_Y(x)分别为与集合X合集合Y同均值同方差的高斯函数；

5b)将集合X和Y的概率密度f_X和f_Y概率密度代入交叉熵方程：

$K\left(Y\right|X)=\∫\log \frac{f_X\left(x\right)}{f_Y\left(x\right)}{f}_X\left(x\right)\mathrm{dx}$

$K\left(X\right|Y)=\∫\log \frac{f_Y\left(x\right)}{f_X\left(x\right)}{f}_Y\left(x\right)\mathrm{dx}---\left(5H\right)$

KLD＝K(Y|X)+K(X|Y)

得到的KLD值作为窗口中心像素点的熵值，此熵值作为变化差异图对应像素点的像素值；

5c)重复步骤5b)，对每一个像素点进行计算，得到变化差异图。

Images