CN102005050B - 用于变化检测的高斯对数模型单边曲率拟合阈值方法 - Google Patents

Abstract

本发明公开了用于变化检测的高斯对数模型单边曲率拟合阈值方法，属于SAR图像变化检测领域。该方法的实现过程为：首先对两幅不同时间相同地域的SAR图像构造差异图并求出差异图像的直方图，在直方图上确定单边拟合区域，假设单边拟合区域中灰度级为阈值，通过单边拟合评价函数计算高斯对数模型对单边拟合区域的曲率拟合误差，求出单边拟合区域中曲率拟合误差最小的阈值，确定无变化区域的直方图概率分布函数和初始阈值，根据初始阈值确定基于高斯模型变化区域直方图概率分布函数。最后用最大后验概率方法确定最终阈值，通过该阈值构造SAR图像变化检测结果图，完成对两幅不同时间、相同地域SAR图像变化区域的检测。通过与几种SAR图像变化检测阈值方法对比，本发明得到的SAR图像变化检测结果的性能最好。

Description

用于变化检测的高斯对数模型单边曲率拟合阈值方法

技术领域

本发明涉及SAR图像变化检测领域，是处理由不同时相SAR图像构造的比值差异影像图的阈值方法。具体地说是一种用于变化检测的高斯对数模型单边曲率拟合阈值方法，来解决SAR图像变化检测领域中变化区域正确检测率不高的问题，本方法适用于多种复杂SAR图像，用于提高SAR图像变化检测系统的检测精度和系统稳定性。

背景技术

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

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

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

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

发明内容

本发明的目的在于：为了克服现有方法具有较高的检测错误率之不足，进一步提高SAR图像变化检测的正确率，从而提高检测精度，针对SAR图像变化检测自身的特点，提出了一种用于变化检测的高斯对数模型单边曲率拟合阈值方法，与其它经典的方法和近来提出的方法相比能获得较高的检测正确率。

实现本发明目的的技术方案是：首先对两幅不同时间、相同地域的SAR图像构造差异图并求出差异图像的直方图，然后在直方图上确定单边拟合区域，接着通过单边拟合评价函数计算高斯对数模型对单边拟合区域的直方图曲率的拟合误差，之后求出单边拟合区域中曲率拟合误差最小的阈值，并确定无变化区域的直方图概率分布函数和初始阈值。然后根据初始阈值确定基于高斯模型的变化区域直方图概率分布函数，最后应用最大后验概率方法确定最终的阈值，并通过该阈值构造SAR图像变化检测结果图，完成对两幅不同时间、相同地域SAR图像变化区域的检测，其具体步骤如下：

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

(2)遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值，求出该差异图像的直方图h(g)，g∈[0，255]为灰度级范围，在直方图上确定单边拟合区域，区域范围的上限为mup，下限为mlow；

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

(4)根据阈值T_i应用高斯对数模型

计算阈值T_i下的无变化区域的直方图概率分布函数p_ncln(x|μ，σ)，式中μ为均值参数，σ为方差参数，x＞0为直方图灰度级；

(5)利用评价函数OSC(g，T_i)＝var(|h(g)-p_ncln(g|μ，σ)|)，g∈[mlow，mup]为直方图灰度级，计算单边曲率拟合误差；

(6)利用公式c(i)＝min(k|p_ncln(x||μ±μ*ln(k)|，σ)，k∈(1，μ])式中μ为均值，k为参考平移量，σ为方差，确定阈值T_i下拟合单边区域曲线的平移量c；

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

(8)求出单边拟合区域中曲率拟合误差最小的阈值T_i，确定无变化区域的直方图概率分布函数和初始阈值；

(9)根据初始阈值，利用高斯模型和最大期望EM参数估计方法，求出变化区域直方图曲线的概率分布函数；

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

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

1、本发明利用SAR图像变化检测特点，在对直方图进行建模时，仅仅对可能产生阈值的区域，即单边区域进行建模，能够更加准确的拟合该区域曲线，提高阈值的精确性。

2、在建立单边区域的模型时，应用高斯对数模型对单边直方图区域的直方图曲线进行曲率拟合，能够很好好地拟合单边直方图曲线。

3、本发明与其它经典的SAR图像变化检测阈值方法相比，能够获得较高的检测正确率，所获得的混合概率分布能够最佳拟合单边区域的直方图分布曲线，验证了本发明的有效性；得到的SAR图像变化检测结果的性能最好。

4、本发明只考虑直方图部分区域的曲线拟合情况，而无需考虑整个直方图曲线的拟合情况，增加了拟合的精确度，从而提高了最终阈值的精度，避免了由于单边区域曲线拟合不精确所导致的阈值偏差，从而使得在SAR图像变化检测中得到更好的结果。

附图说明

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

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

图3是本发明涉及的墨西哥市郊区域两幅SAR影像图、标准图和差异图，其中，(a)表示2002年5月该地区地貌信息，(b)表示2005年4月该地区地貌信息，(c)表示变化检测的标准结果图，(d)表示通过(a)(b)两图SAR图像构造的差异图像

图4是本发明方法和三种对比方法对墨西哥市郊区的实验结果图，其中，(LN)代表基于Log-Normal分布的广义K&I阈值方法的实验结果图，(NR)代表基于Nakagami分布的广义K&I阈值方法的实验结果图，(WR)代表基于Weibull分布的广义K&I阈值方法的实验结果图，(SLNT)代表本发明阈值方法的实验结果图

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

图6是本发明方法和三种对比方法农田区域差异图的单边直方图拟合比较图，其中，(a)代表基于Log-Normal分布的广义K&I阈值方法对农田区域SAR图像数据直方图的拟合情况。(b)代表该方法单边区域直方图的局部放大拟合情况(c)代表基于Nakagami分布的广义K&I阈值方法对农田区域SAR图像直方图拟合情况。(d)代表该方法单边区域直方图的局部放大拟合情况(e)本发明方法对农田区域SAR图像直方图拟合情况。(f)代表本方法单边区域直方图的局部放大拟合情况

图7是本发明方法和三种对比方法对农田区域的实验结果图，其中，(SLNT)代表本发明阈值方法的实验结果图，(NR)代表基于Nakagami分布的广义K&I阈值方法的实验结果图，(LN)代表基于Log-Normal分布的广义K&I阈值方法的实验结果图

具体实施方式

实现本发明一种用于变化检测的高斯对数模型单边曲率拟合阈值方法，首先对两幅不同时间相同地域的SAR图像构造差异图并求出差异图像的直方图，然后在直方图上确定单边拟合区域，接着假设单边拟合区域中的灰度级为阈值，通过单边拟合评价函数计算高斯对数模型对单边拟合区域的曲率拟合误差，之后求出单边拟合区域中曲率拟合误差最小的阈值，并确定无变化区域的直方图概率分布函数和初始阈值。然后根据初始阈值确定基于高斯模型的变化区域直方图概率分布函数。最后应用最大后验概率方法确定最终的阈值，并通过该阈值构造SAR图像变化检测结果图，完成对两幅不同时间、相同地域SAR图像变化区域的检测。下边参照图1，对该发明的实现过程进行详细说明：

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

通过公式

构造比值差异图像，使得在比值差异图像上低灰度级代表无变化区域，高灰度级代表变化区域，然后应用李滤波方法对比值差异图像去除噪声，得到差异图像。

2、遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值，求出该差异图像的灰度直方图h(g)，g∈[0，255]为灰度级范围，在直方图上确定单边拟合区域，区域范围的上限为mup，下限为mlow。

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

(2)参考图2其中，(a)为理想变化检测直方图，(b)为非理想变化检测直方图(无混叠)，(c)为非理想变化检测直方图(混叠)，总结了SAR图像变化检测中差异图的三种直方图分布情况，而最终的阈值都是在Mn和Mc之间产生，称Mn和Mc之间的区域为单边拟合区域，区域范围的上限为mup，下限为mlow。其范围可通过下面公式计算得出。

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

mup＝min(i||h(x)-h(x+1)|＜thd，x∈[mlow，L-2])    (2)

其中，thd＝10^-6，x为直方图灰度级。

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

初始化阈值T_i，从单边区域下界mlow开始进行最佳阈值的选取。

4、根据阈值T_i应用高斯对数模型计算阈值T_i下的无变化区域的直方图概率分布函数p_ncln(x|μ，σ)，式中，x＞0为直方图灰度级，式中μ为均值参数，σ为方差参数，x＞0为直方图灰度级；

应用高斯对数模型

$p\left(x\right|\mathrm{\μ},\mathrm{\σ})=\frac{1}{\mathrm{x\σ}\sqrt{2\mathrm{\π}}}\exp \{-\frac{{(\ln x-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}\},x>0---\left(3\right)$

计算阈值T_i下的无变化区域的直方图概率分布函数p_ncln(x|μ，σ)，μ，σ分别为均值和方差参数。如式(3)可以看出，高斯对数模型分布需要进行估计的参数主要有μ、σ，则通过无变化区域的直方图μ和σ参数估计方法如下：

$\mathrm{\μ}=\frac{\underset{x\∈{\mathrm{\ω}}_n}{\mathrm{\Σ}}h\left(x\right)\ln x}{\underset{x\∈{\mathrm{\ω}}_n}{\mathrm{\Σ}}h\left(x\right)}---\left(4\right)$

$\mathrm{\σ}=\frac{\underset{x\∈{\mathrm{\ω}}_n}{\mathrm{\Σ}}h\left(x\right){(\ln x-\mathrm{\μ})}^2}{\underset{x\∈{\mathrm{\ω}}_n}{\mathrm{\Σ}}h\left(x\right)}---\left(5\right)$

其中，ω_n代表无变化区域，x＞0为直方图灰度级，h(x)为直方图曲线。

5、利用评价函数OSC(g，T_i)＝var(|h(g)-p_ncln(g|μ，σ)|)，g∈[mlow，mup]为直方图灰度级，计算单边曲率拟合误差。

在单边拟合区域内，当该区域的每个灰度级为阈值T_i时，通过步骤4求出相应阈值下的无变化区域的直方图概率分布函数p_ncln(x|μ，σ)，因为在本方法中拟合的是单边曲线，是一种局部直方图的拟合方法，所以本方法中利用p_ncln(x|μ，σ)拟合单边曲线的曲率，故构造单边拟合评价函数OSC(g，T_i)＝var(|h(g)-p_ncln(g|μ，σ)|)，g∈[mlow，mup]计算出单边区域内每个灰度级为阈值时，无变化区域概率分布函数拟合单边区域直方图曲率的情况。

6、利用公式c(i)＝min(k|p_ncln(x||μ±μ*ln(k)|，σ)，k∈(1，μ])，式中μ为均值，k为参考平移量，σ为方差，确定阈值T_i下拟合单边区域曲线的平移量c。

本方法中的单边拟合方法，是通过无变化区域的直方图概率分布函数拟合直方图单边区域，采用的是曲率拟合，所以最佳的分布函数在单边区域内是曲率的最佳拟合，所以对分布函数进行一定的平移才能完全拟合单边区域的曲线，于是本方法构造公式c(i)＝min(k|p_ncln(x||μ±μ*ln(k)|，σ)，k∈(1，μ])μ为均值，利用该公式阈值T_i下拟合单边区域曲线的平移量c。

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

判断是否满足终止条件T_i≤mup，不满足则转步骤4否则转步骤8

8、求出单边拟合区域中曲率拟合误差最小的阈值T_i，确定无变化区域的直方图概率分布函数和初始阈值。

在单边拟合区域中求出拟合误差最小的阈值T_i，通过下列公式

i＝min(k|OSC(g，T_k)，T_k∈[mlow，mup])    (6)

其中，k为单边区域中对应的阈值标号，i为曲率拟合误差最小的阈值对应的标号。同时确定了平移量c(i)，从而确定了初始阈值T_c＝T_i+c(i)，然后通过阈值T_i和对应的平移量c(i)确定无变化区域的直方图概率分布函数。步骤如下：

(1)通过公式(4)(5)分别确定参数μ和σ。

(2)然后将μ-c(i)和σ代入公式(3)中(即p_ncln(x|μ-c(i)，σ))求得无变化区域的直方图概率分布函数。

9、根据初始阈值，利用高斯模型和最大期望(EM)参数估计方法，求出变化区域直方图曲线的概率分布函数。

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

$p_{\mathrm{guass}}\left(x\right|{\mathrm{\ω}}_c)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}-\exp \left(\frac{-{(x-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}\right)---\left(7\right)$

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

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

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

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

${\mathrm{\ω}}_k=\arg \underset{{\mathrm{\ω}}_i\∈\{{\mathrm{\ω}}_n,{\mathrm{\ω}}_c\}}{\max }\left\{p\left({\mathrm{\ω}}_i\right|X(i,j))\right\}=\arg \underset{{\mathrm{\ω}}_i\∈\{{\mathrm{\ω}}_c,{\mathrm{\ω}}_c\}}{\max }\left\{P\left({\mathrm{\ω}}_i\right)p\left(X(i,j)\right|{\mathrm{\ω}}_i)\right\}---\left(8\right)$

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

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

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

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

本发明方法的性能通过如下两个计算机仿真实验给出。

为了验证用于变化检测的高斯对数模型单边曲率拟合阈值方法的优越性，我们将其与基于Nakagami分布、Log-Normal分布和Weibull分布的三种广义K&I阈值算法性能做比较。分别对带参考图的真实SAR图像数据：墨西哥市郊SAR图像数据，以及不带参考图的黄河入海口地区的农田区域SAR图像数据进行了实验。这里把“基于Nakagami分布的广义K&I阈值方法”简称为“NR”，把“基于Log-Normal分布的广义K&I阈值方法”简称为“LN”，把“基于Weibull分布的广义K&I阈值方法”简称为“WR”，把“用于变化检测的高斯对数模型单边曲率拟合阈值方法”简称为“SLNT”。

实验1.带参考图的真实墨西哥市郊区域SAR图像数据

这里分别应用各种算法在真实墨西哥市郊区域SAR图像数据上。实验相关图像如参考图3所示。

(a)表示2002年5月该地区地貌信息，(b)表示2005年4月该地区地貌信息，(c)表示变化检测的标准结果图，(d)表示通过(a)(b)两图SAR图像构造的差异图像。

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

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

为了体现各种算法的性能如下表所示：

表1  墨西哥市郊区域各种算法变化检测性能结果

从表1中可以看出，本发明方法在变化检测总错误数上，比其他对比实验方法都要少，而且在误检数与漏检数上相对都比较均衡。同时从参考图4中可看出，本发明在视觉效果上也较好，其他方法要么噪声较大要么轮廓不清晰，说明都存在欠分割和过分割的现象，体现了本发明方法良好的性能。总体来看，通过表1性能结果和参考图4的视觉效果说明了本发明方法对单边区域的拟合能够达到很高的拟合精度，从而保证了阈值的准确确定，提高了最终变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。其他对比算法，由于是对整体直方图区域的拟合策略，而且依赖于对直方图模型的精确确定，所以一旦整体上较大程度的拟合了直方图曲线，但在单边区域又不能很精确的拟合单边区域直方图曲线，降低了最终阈值准确度，影响到最终的变化检测结果。综上所述本发明方法能够得到更好的SAR图像变化检测结果。

实验2.不带参考图的真实农田区域SAR图像数据

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

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

各种算法在真实农田区域SAR图像数据的单边区域直方图拟合情况如参考图6所示。

在参考图6中，(a)代表基于Log-Normal分布的广义K&I阈值方法对农田区域SAR图像数据直方图的拟合情况。(b)代表该方法单边区域直方图的局部放大拟合情况(c)代表基于Nakagami分布的广义K&I阈值方法对农田区域SAR图像直方图拟合情况。(d)代表该方法单边区域直方图的局部放大拟合情况(e)本发明方法对农田区域SAR图像直方图拟合情况。(f)代表本方法单边区域直方图的局部放大拟合情况。

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

在参考图7中，(SLNT)代表本发明阈值方法的实验结果图，(NR)代表基于Nakagami分布的广义K&I阈值方法的实验结果图，(LN)代表基于Log-Normal分布的广义K&I阈值方法的实验结果图。

从参考图7中可以看出，本发明方法轮廓很清新，而且含有噪声也较少，其对变化区域的检测最为准确。同时与LN方法的变化检测结果相比较，LN方法检测出来的变化区域较模糊，很多变化的像素点没有检测出来。与NR方法的变化检测结果相比较，很明显NR方法检测结果的噪声较大，变化检测的结果较差，产生以上结果的原因可以从参考图6中得到答案，因为其他对比方法的单边区域并没有很好的拟合，使阈值产生了偏差。而本发明方法很好地拟合了单边区域的直方图曲线，所以变化检测结果最理想。在本实验中由于WR方法并没有给出合理的结果，所以没有该方法的实验结果。同时体现了WR方法对某些SAR图像不适合，说明了WR方法的局限性。综合参考图6、7及以上分析，可以得出本发明方法的结果与其他方法相比更为精确，单边区域拟合的好坏直接影响着分割阈值的结果，从而直接影响SAR图像变化检测的结果。

本发明提出的方法将高斯对数模型应用到产生阈值的单边区域，进行单边曲线的曲率拟合，并没有像其他方法一样，进行总体的直方图拟合，而是进行局部直方图的曲率拟合，大大增加了拟合的精度，从而增加了SAR图像变化检测的检测精度。

Claims (1)

1.用于变化检测的高斯对数模型单边曲率拟合阈值方法，其特征在于：首先对两幅不同时间、相同地域的SAR图像构造差异图并求出差异图像的直方图，在直方图上确定单边拟合区域，假设单边拟合区域中的灰度级为阈值，通过单边拟合评价函数计算高斯对数模型对单边拟合区域的曲率拟合误差，求出单边拟合区域中曲率拟合误差最小的阈值，并确定无变化区域的直方图概率分布函数和初始阈值，根据初始阈值确定基于高斯模型的变化区域直方图概率分布函数，应用最大后验概率方法确定最终的阈值，并通过该阈值构造SAR图像变化检测结果图，完成SAR图像变化区域的检测，其具体步骤如下：

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

(2)遍历差异图像的每个像素，记录图像上每个灰度级出现的频率值，求出该差异图像的灰度直方图h(g)，g∈[0，255]为灰度级范围，在直方图上确定单边拟合区域，区域范围的上限为mup，下限为mlow；

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

(4)根据阈值T_i应用高斯对数模型 $p\left(x\right|\mathrm{\μ},\mathrm{\σ})=\frac{1}{\mathrm{x\σ}\sqrt{2\mathrm{\π}}}\exp \{-\frac{{(\ln x-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}\}$ 计算阈值T_i下的无变化区域的直方图概率分布函数p_ncln(x|μ，σ)，式中μ为均值参数，σ为方差参数，x＞0为直方图灰度级；

(5)利用评价函数OSC(g，T_i)＝var(|h(g)-p_ncln(g|μ，σ)|)，g∈[mlow，mup]为直方图灰度级，计算单边曲率拟合误差；

(6)利用公式c(i)＝min(k|p_ncln(x||μ±μ*ln(k)|，σ)，k∈(1，μ])，式中μ为均值，k为参考平移量，σ为方差，确定阈值T_i下拟合单边区域曲线的平移量c；

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

(8)求出单边拟合区域中曲率拟合误差最小的阈值T_i，并确定无变化区域的直方图概率分布函数和初始阈值；

(9)根据初始阈值，利用高斯模型和最大期望EM参数估计方法，求出变化区域直方图曲线的概率分布函数；

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

Images