CN102005049B - 基于单边广义高斯模型的sar图像变化检测阈值方法 - Google Patents

Abstract

本发明公开一种基于单边广义高斯模型的SAR图像变化检测阈值方法，属于SAR图像变化检测领域。其实现过程为：首先对两幅不同时间相同地域的SAR图像构造差异图，求出差异图像的直方图，接着利用单边广义高斯模型求出其无变化区域直方图概率分布函数，利用高斯模型求出其变化区域直方图概率分布函数，通过最大后验概率方法自动确定阈值，再通过该阈值生成变化检测结果图，最终检测出两幅SAR图像中变化区域。本发明利用新构造模型对阈值产生区域直方图曲线进行拟合，提高了最终阈值精度，避免了曲线拟合不精确导致的阈值偏差，使SAR图像变化检测结果的性能最佳。

Description

基于单边广义高斯模型的SAR图像变化检测阈值方法

技术领域

本发明属于SAR图像变化检测领域，涉及SAR图像变化检测中的阈值技术。具体地说是提出了一种基于单边广义高斯模型的阈值方法，用来解决SAR图像变化检测领域中变化区域检测错误率较高的问题，提高SAR图像变化检测中的检测精度及速度。

背景技术

SAR图像变化检测是从不同时间获取同一地理区域的多时相遥感影像，定性地分析和确定地表变化过程和特征的技术。由于与光学遥感系统相比，SAR系统具有全天时、全天候获取数据的能力，所以SAR图像变化检测技术正广泛的应用于各个领域，例如环境监控，农业研究，城市区域研究，森林监控等方面。

阈值技术是SAR图像变化检测中的关键技术之一。该技术通过变化前后不同时相的SAR影像图构造差异图，然后在差异图的基础上进行自动阈值的确定，形成变化检测结果图，完成变化检测过程，该技术具有方法简单，时间迅速等特点。同时，在相关技术中该技术又往往作为其中的某个环节，其精度高低也直接影响相关算法的性能。因此，该技术在SAR图像变化检测领域中至关重要，得到了广泛的关注。在阈值技术中，差异影像直方图的概率统计分布是其核心内容，概率统计模型能否很好的拟合差异影像直方图，直接影响着变化检测错误率的大小，受到了国内外研究者的广泛关注。所以，目前的研究主要集中在统计概率分布模型的建立上。

国内外学者提出了不同的阈值方法来提高SAR图像变化检测系统的检测性能，但均存在不同的问题，主要检测错误率较高的问题。阈值方法作为一种实时性较好方法，具有自动性，实时性的优点，但是现有的阈值方法检测错误率仍较高。

在SAR图像变化检测中，许多学者给出了许多不同的阈值方法。Kittler和Illingworth提出了经典了最小误差阈值方法，该方法假设变化区域与非变化区域的直方图统计分布符合高斯分布，通过最小化惩罚函数来自动确定阈值。但该方法中基于高斯模型的假设并不准确，实际SAR图像数据分布并不符合高斯分布，所以其检测错误率仍较高。意大利G.Moser，S.B.Serpico等人在最小误差阈值基础上提出了一种广义最小误差阈值方法，该方法以三种概率分布模型：Nakagami分布、Log-Normal分布和Weibull分布为基础扩展了最小误差阈值方法。该阈值方法由于使用了更加适合SAR图像数据分布的模型，所以呈现了较好的检测结果。但是，这三种模型仍不能完全拟合差异图直方图，所以其检测错误率仍然较高。

发明内容

本发明的目的在于：为了降低SAR图像变化检测的错误率克服现有方法较高的检测错误率，提高检测精度，针对SAR图像变化检测的特点，提出了基于单边广义高斯模型的阈值方法，与其它现有的方法相比能够获得较低的检测错误率。

本发明的技术方案是：首先对两幅不同时间相同地域的SAR图像构造差异图，然后求出差异图像的直方图，接着利用单边广义高斯模型求出其无变化区域的直方图概率分布函数，再利用高斯模型求出其变化区域的直方图概率分布函数，最后通过最大后验概率方法自动确定阈值，再通过该阈值生成变化检测结果图，完成对两幅SAR图像中变化区域的最终检测，其具体实现步骤如下：

(1)对两幅不同时间相同地域的SAR图像I₁，I₂构造差异图像；

(2)遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值，求出该差异图像的直方图h(z)，z∈[0，255]为灰度级范围，在直方图上确定可能产生阈值的灰度级区域，即阈值产生区域，该区域的上限为m_n_up，下限为m_n_low；

(3)令阈值T_i＝m_n_low，i＝1；

(4)构造单边广义高斯模型p(X)＝aexp{-[b|X-(m_n-c)|]^β}构造该模型相对应的参数估计方法，求出阈值T_i下的无变化区域直方图曲线的概率分布函数p_nosgg(z|T_i)，其中参数σ是方差参数，β是形状控制参数，m_n是均值参数，c是单边平移参数，X为直方图灰度级，Γ(·)为Gamma函数，其表达式为

z为直方图灰度级；

(5)利用准则函数

计算在阈值T_i下无变化区域直方

图概率分布函数p_nosgg(z|T_i)所对应的曲线拟合误差；

(6)令T_i+1＝T_i+1，i＝i+1，如果T_i+1≤m_n_up，转步骤(4)，否则转步骤(7)；

(7)令阈值产生区域中拟合误差最小的阈值为初始阈值T_c，同时确定无变化区域直方图曲线的概率分布函数；

(8)利用高斯模型与期望最大化EM参数估计方法，求出变化区域直方图曲线的概率分布函数；

(9)根据最大后验概率方法确定最终阈值，并通过该阈值构造变化检测结果图。

本发明与现有的技术相比具有如下有益效果：

1、本发明在广义高斯模型的基础上构造了一种单边广义高斯模型，能够更好地拟合直方图阈值产生区域的曲线。从而提高了最终阈值的精度，避免了由于阈值产生区域曲线拟合不精确所导致的阈值偏差，从而使得SAR图像变化检测得到更好的结果。通过与几种SAR图像变化检测阈值方法对比，本发明方法得到的SAR图像变化检测结果的性能最佳。

2、本发明在应用单边广义高斯模型拟合曲线时，只考虑阈值产生区域直方图曲线的拟合情况，而无需考虑整个直方图的拟合情况，增加了阈值产生区域曲线的拟合精确度，从而保证了阈值的精确确定。

3、本发明与其它经典的SAR图像变化检测阈值方法相比，获得了较低的检测错误率，所获得的混合概率分布曲线在单边直方图区域能够达到最佳的拟合程度，验证了本发明的有效性。

附图说明

图1是本发明方法的流程图

图2是本发明概括的三种不同特性SAR图像变化检测直方图，其中，(a)为理想变化检测直方图，为非理想变化检测直方图(无混叠)，(c)为非理想变化检测直方图(混叠)

图3是本发明涉及的Ottawa地区两幅SAR影像图、标准图和差异图，(a)表示1997年5月Ottawa地区的地貌信息。(b)表示1997年8月Ottawa地区的地貌信息。(c)表示变化检测的标准结果图(d)表示通过(a)(b)两图SAR图像构造的差异图像

图4是本发明方法和三种对比方法对Ottawa地区SAR影像的实验结果图，其中，(LN)代表基于Log-Normal分布的广义最小误差阈值方法的实验结果图，(NR)代表基于Nakagami分布的广义最小误差阈值方法的实验结果图，(WR)代表基于Weibull分布的广义最小误差阈值方法的实验结果图，(SGGT)代表本发明方法的实验结果图

图5本发明两幅不同时相黄河入海口区域两幅SAR影像图和差异图，其中，(a)表示2008年6月入海口区域地貌信息，(b)表示2009年6月入海口区域地貌信息，(c)表示通过(a)(b)两图SAR图像构造的差异图像

图6是本发明方法和三种对比方法入海口区域比值差异图的直方图拟合比较图，其中(a)表示本发明方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况，(b)表示LN-GKIT方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况，(c)表示NR-GKIT方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况，(d)表示WR-GKIT方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况

图7是本发明方法和三种对比方法对入海口区域SAR影像的实验结果图，其中(SGGT)代表本发明阈值方法的实验结果图，(NR-GKIT)代表基于Nakagami分布的广义最小误差阈值方法的实验结果图，(WR-GKIT)代表基于Weibull分布的广义最小误差阈值方法的实验结果图，(LN-GKIT)代表基于Log-Normal分布的广义最小误差阈值方法的实验结果图

具体实施方式

参照图1，实现本发明一种基于单边广义高斯模型的SAR图像变化检测阈值方法，首先对两幅不同时间相同地域的SAR图像构造差异图，然后求出差异图像的直方图，接着利用单边广义高斯模型求出其无变化区域的直方图概率分布函数，再利用高斯模型求出其变化区域的直方图概率分布函数，最后通过最大后验概率方法自动确定阈值，再通过该阈值生成变化检测结果图。下边对该发明的实现过程进行详细说明：

1、对两幅不同时间相同地域的SAR图像I₁，I₂构造差异图像：

构造方法通过公式

构造比值差异图像DI，使得在比值差异图像上低灰度级呈现为无变化区域，高灰度级呈现为变化区域，然后应用均值滤波对比值差异图像去除噪声，得到最终差异图像；

2、遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值，求出该差异图像的直方图h(z)，z∈[0，255]为灰度级范围，在直方图上确定可能产生阈值的灰度级区域，即阈值产生区域，该区域的上限为m_n_up，下限为m_n_low：

(1)遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值求得差异图像的灰度直方图h(z)，z∈[0，255]为灰度级范围；

(2)参考图2总结了SAR图像变化检测中差异图的三种直方图分布情况，而最终的阈值都是在Mn和Mc之间产生，称Mn和Mc之间的区域为阈值产生区域，区域的范围为[m_n_low，m_n_up]。其范围可通过下面公式计算得出：

m_n_low＝max(h(x))x∈[0，255]             (1)

$m_{n\_}\mathrm{up}=\min \left(i\right|\left|\frac{\underset{C=1}{\overset{T}{\mathrm{\&Sigma;}}}h(i+c)-\underset{C=2}{\overset{T+1}{\mathrm{\&Sigma;}}}h(i+c)}{T}\right|<e,i\&Element;[m_{n\_}\mathrm{low},L-1-T])---\left(2\right)$

其中，i为直方图灰度级，c为平滑因子步长，e＝10^-6，T为直方图平滑因子，取T＝10；

3、令阈值T_i＝m_n_low，i＝1；

初始化阈值T_i＝m_n_low，i＝1，在阈值产生区域寻找最佳阈值；

4、构造单边广义高斯模型p(X)＝aexp{-[b|X-(m_n-c)|]^β}构造该模型相对应的参数估计方法，求出阈值T_i下的无变化区域直方图曲线的概率分布函数p_nosgg(z|T_i)，其中

参数σ是方差参数，β是形状控制参数，m_n是均值参数，c是单边平移参数，X为直方图灰度级，Γ(·)为Gamma函数，其表达式为

z为直方图灰度级；

构造单边广义高斯模型如下：

p(X|ω_n)＝aexp{-[bX-(m_n-c)|^β}X∈[0，T_i]         (3)

$\left\{\begin{array}{c}a=\frac{\mathrm{b\&beta;}}{2\mathrm{\&Gamma;}(1/\mathrm{\&beta;})}\\ b=\frac{1}{\mathrm{\&sigma;}}\sqrt{\frac{\mathrm{\&Gamma;}(3/\mathrm{\&beta;})}{\mathrm{\&Gamma;}(1/\mathrm{\&beta;})}}\\ m_n=\frac{\underset{X=0}{\overset{T_i}{\mathrm{\&Sigma;}}}X\&CenterDot;h\left(X\right)}{\underset{X=0}{\overset{T_i}{\mathrm{\&Sigma;}}}h\left(X\right)}\end{array}\right.---\left(4\right)$

式(3)和(4)中，参数σ，β，m_n，c分别是单边广义高斯分布的方差、形状控制、均值和单边平移参数。式(4)中的Γ(·)为Gamma函数，其表达式为

假设阈值产生区域内的灰度级都为阈值，计算阈值T_i下的无变化区域的直方图概率分布函数p_nosgg(z|T_i)；

从式(3)可以看出，单边广义高斯模型需要进行估计的参数主要有c，β，则通过无变化区域的直方图，β和c参数估计方法如下：

(1)参数β估计步骤如下：

[1]通过计算所谓的广义高斯率函数r(β)产生一个查找表，其中r(β)如下式：

$r\left(\mathrm{\&beta;}\right)=\frac{\mathrm{\&Gamma;}\left(\frac{1}{\mathrm{\&beta;}}\right)\mathrm{\&Gamma;}\left(\frac{3}{\mathrm{\&beta;}}\right)}{{\mathrm{\&Gamma;}}^2\left(\frac{2}{\mathrm{\&beta;}}\right)}---\left(5\right)$

[2]估计均值的绝对期望E[z|ω_n]。

$E\left[z\right|{\mathrm{\&omega;}}_n]=\frac{1}{\underset{X=0}{\overset{T_i}{\mathrm{\&Sigma;}}}h\left(X\right)}\underset{z=0}{\overset{T_i}{\mathrm{\&Sigma;}}}h\left(z\right)\&CenterDot;|z-m_n|---\left(6\right)$

其中，ω_n代表无变化区域，参数z，m_n分别为直方图灰度级和均值；

[3]计算比率ρ_n

${\mathrm{\&rho;}}_n=\frac{{\mathrm{\&sigma;}}_n^2}{E^2\left[z\right|{\mathrm{\&omega;}}_n]}---\left(7\right)$

其中，σ_n为方差参数；

[4]通过查找表计算

等式的解。利用反比例函数最小二乘拟合求得r^-1(ρ_u)，即

$r^{-1}\left(x\right)=\frac{0.2718}{0.7697-x}-0.1247---\left(8\right)$

(2)参数c估计步骤如下：

[1]初始化c＝0，将(1)中估计的参数β代入式(4)得到a，b，将参数βabc代入式(3)中；

[2]令

End＝m_n_up，则

$c=\frac{\underset{z=\mathrm{Start}}{\overset{\mathrm{End}}{\mathrm{\&Sigma;}}}{p}^{-1}\left(h\left(z\right)\right|{\mathrm{\&omega;}}_n)-z}{\mathrm{Start}-\mathrm{End}+1}---\left(9\right)$

其中，p^-1(x|ω_n)为阈值产生区域曲线所对应的广义高斯模型的反函数，可以通过建立查找表来求解该方程的解；

5、利用准则函数

z为直方图灰度级，计算在阈值T_i下无变化区域直方图概率分布函数p_nossgg(z|T_i)所对应的曲线拟合误差；

在阈值产生区域内，当该区域的每个灰度级为阈值T_i时，通过步骤4求出相应阈值T_i下的无变化区域直方图概率分布函数p_nosgg(z|T_i)，然后利用准则函数

计算在阈值T_i下无变化区域直方图概率分布函数p_nosgg(z|T_i)所对应的曲线拟合误差；

6、令T_i+1＝T_i+1，i＝i+1，如果T_i+1≤m_n_up，转步骤4，否则转步骤7；

判断是否满足终止条件T_i+1≤m_n_up，不满足则转步骤4否则转步骤7；

7、令阈值产生区域中拟合误差最小的阈值为初始阈值T_c，同时确定无变化区域直方图曲线的概率分布函数；

通过准则函数C(z，T_i)在阈值产生区域找到拟合误差最小的阈值，并令该阈值为初始阈值T_c，初始阈值T_c确定后，通过步骤4求出无变化区域直方图曲线的概率分布函数；

8、利用高斯模型与期望最大化EM参数估计方法，求出变化区域直方图曲线的概率分布函数；

假设整个变化区域的概率分布函数符合高斯分布p(X|ω_c)如下：

$p_{\mathrm{guass}}\left(X\right|{\mathrm{\&omega;}}_c)=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \left(\frac{-{(x-\mathrm{\&mu;})}^2}{2{\mathrm{\&sigma;}}^2}\right)---\left(10\right)$

在式(10)中，参数σ与μ分别是方差和均值，ω_c代表变化区域直方图灰度级；

通过最大期望EM参数估计方法估计出参数σ与μ，则确定了变化区域直方图曲线的概率分布函数；

9、根据最大后验概率方法确定最终阈值，并通过该阈值构造变化检测结果图；

根据最大后验概率方法，在差异图中，每个像素点X(i，j)属于哪类可根据下式给出；

${\mathrm{\&omega;}}_k=\arg \underset{{\mathrm{\&omega;}}_i\&Element;\{{\mathrm{\&omega;}}_n,{\mathrm{\&omega;}}_c\}}{\max }\left\{p\left({\mathrm{\&omega;}}_i\right|X(i,j))\right\}=\arg \underset{{\mathrm{\&omega;}}_i\&Element;\{{\mathrm{\&omega;}}_n,{\mathrm{\&omega;}}_c\}}{\max }\left\{P\left({\mathrm{\&omega;}}_i\right)p\left(X(i,j)\right|{\mathrm{\&omega;}}_i)\right\}---\left(11\right)$

应用该准则解决变化检测问题相当于在变化区域和无变化区域求差异图的阈值，因此最终阈值可通过下式得出。

P(ω_c)p(X|ω_c)＝P(ω_n)p(X|ω_n)              (12)

其中P(ω_n)，P(ω_c)为变化区域与无变化区域的先验概率。

根据最终阈值遍历整个图像，像素小于阈值的置0，像素大于等于阈值的置255，则构造出变化检测结果图；

本发明将应用在SAR图像变化系统中，其性能可通过下面两个计算机仿真实验给出。

为了验证基于单边广义高斯阈值方法的优越性，我们将本发明方法与近年来提出的基于Nakagami分布、基于Log-Normal分布和基于Weibull分布的三种广义最小误差阈值算法性能做出比较。分别对带参考图的真实SAR图像数据：Ottawa地区SAR图像数据，以及不带参考图的黄河入海口地区的部分区域SAR图像数据进行了实验。这里把基于Nakagami分布的广义最小误差阈值方法简称为“NR”或“NR-GKIT”，把基于Log-Normal分布的广义最小误差阈值方法简称为“LN”或“LN-GKIT”，把基于Weibull分布的广义最小误差阈值方法简称为“WR”或“WR-GKIT”，把基于单边广义高斯模型SAR图像变化检测阈值方法简称为“SGGT”。

实验1.带参考图的真实Ottawa SAR图像数据

这里分别应用各种方法在真实Ottawa SAR图像数据上。实验相关图像如参考图3所示。(a)表示1997年5月Ottawa地区的地貌信息。(b)表示1997年8月Ottawa地区的地貌信息。(c)表示变化检测的标准结果图(d)表示通过(a)(b)两图SAR图像构造的差异图像。

各种算法在真实Ottawa区域SAR图像数据的实验结果图如参考图4所示。

在参考图4中(LN)代表基于Log-Normal分布的广义最小误差阈值方法的实验结果图，(NR)代表基于Nakagami分布的广义最小误差阈值方法的实验结果图，(WR)代表基于Weibull分布的广义最小误差阈值方法的实验结果图，(SGGT)代表本发明方法的实验结果图。

各种方法的性能指标如下表所示：

表1Ottawa地区各种算法变化检测结果

从表1中可以看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，而且在误检数与漏检数上相对都比较均衡，数量上也处于中间水平，说明了本发明方法的优越性。从参考图4中可以略微看出，LN方法与NR方法的噪声点较多，从表1中的误检数也可以给出同样的结论，这两种方法的误检数是四种方法中最高的。而WR方法有一些变化细节没有体现出来，从表1中的漏检测数也可以给出同样的结论，WR方法的漏检测数是四种方法中最高的。总体来看，通过表1性能结果和参考图4的视觉效果说明了本发明方法对阈值产生区域的拟合能够达到很高的拟合精度，从而保证了阈值的确定，提高了最终变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。相比其他对比方法，由于是对整体直方图区域的拟合策略，所以一旦整体上较大程度的拟合了直方图曲线，但在阈值产生区域又不能很精确的拟合直方图曲线，降低了最终阈值准确度，影响到最终的变化检测结果。综上所述本发明方法能够得到更好的SAR图像变化检测结果。

实验2.不带参考图的真实入海口区域SAR图像数据

这里分别应用各种方法在真实入海口区域SAR图像数据上。实验相关图像如参考图5所示。

(a)表示2008年6月入海口区域地貌信息，(b)表示2009年6月入海口区域地貌信息，(c)表示通过(a)(b)两图SAR图像构造的差异图像。

各种算法在真实入海口区域SAR图像数据的直方图拟合情况和阈值产生区域直方图拟合情况，如参考图6所示。

(a)表示本发明方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况，(b)表示LN-GKIT方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况，(c)表示NR-GKIT方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况，(d)表示WR-GKIT方法对差异图的整体直方图拟合情况和阈值产生区域直方图的拟合情况。

各种算法在真实入海口区域SAR图像数据的实验结果图如参考图7所示。

(SGGT)代表本发明阈值方法的实验结果图，(NR-GKIT)代表基于Nakagami分布的广义最小误差阈值方法的实验结果图，(WR-GKIT)代表基于Weibull分布的广义最小误差阈值方法的实验结果图，(LN-GKIT)代表基于Log-Normal分布的广义最小误差阈值方法的实验结果图。

从参考图7中可以看出，本发明方法虽然含有一定的噪声，但其对变化区域的检测更为准确。与LN-GKIT方法相比变化检测结果相差不多，但LN-GKIT方法变化的像素点稍微少些，漏检测数稍多些。很明显NR-GKIT方法的噪声最大，误检测数最多，变化检测的结果也较差。产生该结果的原因可以从参考图6中得到答案，因为阈值产生区域的直方图并没有很好的拟合，使得阈值产生了偏差。而WR-GKIT方法明显很多变化的区域没有检测出来，主要也是因为单边区域没有较好的拟合，从而导致阈值产生偏差。从参考图6中可以看出，只有本发明方法对阈值产生区域的直方图得到了很好的拟合，而其他方法的拟合情况均较差，主要也是由于其方法整体拟合的策略和其方法自身的缺陷。综合参考图6、7分析，可以得出本发明方法的结果比LN-GKIT方法、NR-GKIT方法和WR-GKIT方法的结果更为精确，单边区域拟合的好坏直接影响着分割的结果，证明了本发明的有效性。

通过两个实验的分析可以看出，本发明方法的结果与其他方法相比更为精确，能够更好地拟合阈值产生区域的直方图，从而得到更佳的SAR图像变化检测的结果。

Claims (1)

1.一种基于单边广义高斯模型的SAR图像变化检测阈值方法，其特征在于：首先对两幅不同时间、相同地域的SAR图像构造差异图，然后求出差异图像的直方图，接着利用单边广义高斯模型求出其无变化区域的直方图概率分布函数，再利用高斯模型求出其变化区域的直方图概率分布函数，最后通过最大后验概率方法自动确定阈值，再通过该阈值生成变化检测结果图，完成对两幅SAR图像中变化区域的最终检测，其具体实现步骤如下：

(1)对两幅不同时间相同地域的SAR图像I₁，I₂构造差异图像；

遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值，求出该差异图像的直方图h(z)，z∈[0，255]为灰度级范围，在直方图上确定可能产生阈值的灰度级区域为m_n_low和m_n_up之间的区域，则m_n_low和m_n_up之间的区域为阈值产生区域，该区域的上限为m_n_up，下限为m_n_low；

(2)令阈值T_i＝m_n_low，i＝1；

(3)构造单边广义高斯模型p(x)＝aexp{-[b|X-(m_n-c)|]^β}构造该模型相对应的参数估计方法，求出阈值T_i下的无变化区域直方图曲线的概率分布函数p_nosgg(z|T_i)，其中

参数σ是方差参数，β是形状控制参数，m_n是均值参数，c是单边平移参数，X为直方图灰度级，Γ(·)为Gamma函数，其表达式为 $\mathrm{\&Gamma;}\left(Z\right)=\underset{0}{\overset{\&infin;}{\&Integral;}}{e}^{-t}{t}^{z-1}\mathrm{dt},$ z为直方图灰度级；

(4)利用准则函数计算在阈值T_i下无变化区域直方图概率分布函数p_nosgg(z|T_i)所对应的曲线拟合误差；

(5)令T_i+1＝T_i+1，i＝i+1，如果T_i+1≤m_n_up，转步骤(4)，否则转步骤(7)；

(6)令阈值产生区域中拟合误差最小的阈值为初始阈值T_c，同时确定无变化区域直方图曲线的概率分布函数；

(7)利用高斯模型与期望最大化EM参数估计方法，求出变化区域直方图曲线的概率分布函数；

(8)根据最大后验概率方法确定最终阈值，并通过该阈值构造变化检测结果图。

Images