CN104778717A - 基于导向差异图的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于导向差异图的SAR图像变换检测方法，主要解决现有技术对同一区域不同时刻的SAR图像之间的变化信息不能精确地描述，及抗干扰能力弱的问题。其实现步骤为：(1)对同一区域不同时刻的两幅SAR图像进行预处理；(2)对预处理后的两幅SAR图像分别构造均值比值差异图和对数比值差异图；(3)根据导向滤波器的思想，将均值比值差异图作为导向图，对数比值差异图作为输入图，进行融合生成导向差异图；(4)利用基于局部相关熵的均值聚类水平集模型对导向差异图进行分割，得到变化检测的结果。本发明能有效地抑制背景区域信息的干扰，提高了变化检测的精确度，可用于自然灾害评估、环境资源监测和城市建设规划。

Description

基于导向差异图的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，具体是涉及一种SAR图像变化检测方法，可用于自然灾害评估、环境资源监测和城市建设规划。

背景技术

合成孔径雷达SAR与可见光和红外遥感相比，具有全天时、全天候、穿透能力强的特性，因此无论在军事领域还是民用领域都得到了广泛的应用。随着SAR技术的快速发展，SAR图像变化检测技术变得日益重要。SAR图像变化检测是指通过比较和分析在同一区域所获得的不同时刻的两幅或多幅SAR图像，根据这些SAR图像之间的变化差异来获取所研究区域随时间而发生的变化信息。SAR图像变化检测的关键步骤在于生成相应的差异图像和对差异图像进一步的分析并提取变化信息。

西安电子科技大学在其专利申请“基于直方图和精英遗传聚类算法的SAR图像变化检测”(专利申请号：CN201310733079.4，公开号：CN103971362A)中提出了一种基于直方图和精英遗传聚类算法的SAR图像变化检测。通过双边滤波器的思想对所生成邻域差值差异图和邻域比值差异图进行融合，能够得到稳定的差异图；但是双边滤波器不具有良好的边缘保持性，而且与对数比差异图像相比，领域差值差异图和邻域比值差异图不能有效地抑制SAR图像的相干斑噪声，造成生成的差异图中变化区域含有较多的干扰信息。

西安电子科技大学在其专利申请“基于改进C-V模型的SAR图像变化检测方法”(专利申请号：CN201210231787.3，公开号：CN102779346A)中提出了一种基于改进C-V模型的SAR图像变化检测方法。通过改进的C-V水平集对差异图进行分割，具有较高的分割精确度，能够有效地提取变化信息；但改进的C-V水平集模型对初始化位置比较敏感且不能有效地抑制噪声的干扰，初始化位置的不同对分割精确度有很大影响，同时具有较高的时间复杂度。

发明内容

本发明的目的在于针对上述现有SAR图像变化检测技术的不足，提出一种基于导向差异图的SAR图像变化检测方法，以有效地抑制差异图像中的干扰信息，增强图像的变化信息，生成鲁棒的差异图像，提高SAR图像变化信息结果提取的精确度。

本发明的技术方案是：通过合成导向差异图和进行LCK水平集进行差异图的分割，提出精确的差异图像分类策略，其实现步骤包括如下：

(1)输入同一区域t₁时刻的SAR图像X₁和t₂时刻的SAR图像X₂进行图像配准、几何校正和图像去噪预处理，获得预处理后两幅不同时刻的SAR图像I₁，I₂；

(2)对预处理后的两幅SAR图像I₁,I₂，通过均值比值法得到均值比值差异图I_M；

(3)对预处理后的两幅SAR图像I₁,I₂，通过对数比值法构造对数比值差异图I_L；

(4)利用导向滤波器，将均值比值差异图I_M作为导向图，将对数比值差异图I_L作为输入图，并对该导向图和输入图进行融合生成导向差异图I_G；

(5)利用基于局部相关熵的K均值聚类水平集模型对导向差异图I_G进行分割，得到变化检测的结果；

(5a)初始化水平集函数为Φ＝Φ⁰，并根据水平集函数Φ在导向差异图I_G上做一条闭合曲线：

${\mathrm{\Φ}}^0=\mathrm{\Φ}(m,t=0)=\left\{\begin{array}{cc}-c& m\∈{\mathrm{\Ω}}_1\\ 0& m\∈{\mathrm{\Ω}}_2\\ c& m\∈{\mathrm{\Ω}}_3\end{array}\right.$

其中，m表示导向差异图I_G中像素点，c表示大于0的常数，c取常数1，Ω₁表示水平集函数Φ将导向差异图I_G分成的闭合曲线内部区域像素点，Ω₂表示水平集函数Φ将导向差异图I_G分成的闭合曲线上像素点，Ω₃表示水平集函数Φ将导向差异图I_G分成的闭合曲线外部区域像素点；

(5b)初始化权重因子：其中，N表示导向差异图I_G的总像素点数，j表示导向差异图I_G中第j个像素点；

(5c)分别计算当前水平集函数Φ将导向差异图I_G分成的内部区域的灰度均值和外部区域的灰度均值t表示表示迭代次数，初始值为0；

(5d)根据内部区域的灰度均值和外部区域的灰度均值通过权重因子的更新公式，更新权重因子

(5e)通过判别规则，对更新后的权重因子进行收缩；

(5f)根据水平集函数的更新公式，对当前水平集函数Φ进行更新；

(5g)对步骤(5c)～(5f)进行重复迭代操作，并判断当前水平集函数Φ是否收敛即Φ^t＝Φ^t-1，如果是，则停止对水平集函数Φ的更新，得到变化检测结果图，完成SAR图像的变化检测；如果否，则继续通过循环迭代对水平集函数进行更新。

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

第一，本发明根据导向滤波器的思想将均值比值差异图作为导向图，对数比值差异图作为输入图进行融合，得到性能良好的导向差异图，使得对同一区域不同时刻的SAR图像之间的变化信息有了更精确地描述。

第二，由于本发明采用LCK水平集模型对导向差异图像中变化区域信息和未变化区域信息进行了精确的分类，有效地避免了噪声信息的干扰，使得变化信息提取的结果更加精确，因此整体提升了SAR图像变化检测的精确度。

附图说明

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

图2是仿真使用的Bern地区不同时刻的两幅SAR图像和变化检测标准图；

图3是用本发明对图2不同时刻的SAR图像得的差异图像；

图4是用本发明与现有方法对Bern地区做的变化检测对比图；

图5是仿真使用的Ottawa地区不同时刻的两幅SAR图像和变化检测标准图；

图6是本发明对图5不同时刻的SAR图像得的差异图像；

图7是本发明与现有方法对Ottawa地区做的变化检测的对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

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

步骤1：对输入图像X₁，X₂预处理。

输入同一区域t₁时刻的SAR图像X₁和t₂时刻的SAR图像X₂，并依次进行图像配准、几何校正和图像去噪预处理，获得预处理后两幅不同时刻的SAR图像I₁，I₂。

步骤2：根据预处理后的图像获得均值比值差异图I_M。

对预处理后的两幅SAR图像I₁，I₂，通过均值比值法构造均值比值差异图I_M，I_M中的每个元素I_M(i,j)按如下公式计算：

$L_M(i,j)=1-\min (\frac{{\mathrm{\ω}}_1(i,j)}{{\mathrm{\ω}}_2(i,j)},\frac{{\mathrm{\ω}}_2(i,j)}{{\mathrm{\ω}}_1(i,j)}),$

其中，ω₁(i,j)表示预处理后t₁时刻的SAR图像I₁在以坐标点为(i,j)处的像素点为中心、邻域窗口大小为5×5的局部均值，ω₂(i,j)表示预处理后t₂时刻的SAR图像I₂在以坐标点为(i,j)处的像素点为中心、邻域窗口大小为5×5的局部均值。

步骤3：根据预处理后的图像获得对数比值差异图I_L。

对预处理后的两幅SAR图像I₁，I₂，通过对数比值法构造对数比值差异图I_L：

$I_L=\left|\log \frac{I_2+1}{I_1+1}\right|=\log (I_2+1)-\log (I_1+1)|.$

步骤4：I_M和I_L融合生成导向差异图I_G。

根据导向滤波器思想，将均值比值差异图I_M作为导向图，将对数比值差异图I_L作为输入图，对这两者进行融合生成导向差异图I_G：

$I_G=\frac{1}{\left|\mathrm{\ω}\right|}\underset{k\backslash i\∈{\mathrm{\ω}}_k}{\mathrm{\Σ}}(a_k{I}_M+b_k)$

其中，|ω|表示ω_k窗口中含有的像素个数，k表示导向差异图I_G的第k个像素点，ω_k表示以第k个像素点为中心、大小为25×25的邻域窗口，i为ω_k邻域窗口中第i个像素点，a_k和b_k为数值不同的线性系数，具体表示如下：

$a_k=\frac{\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}{I}_M\left(i\right)I_L\left(i\right)-{\mathrm{\μ}}_k\stackrel{\‾}{I_L\left(k\right)}}{{\mathrm{\σ}}_k^2+\mathrm{\ϵ}}$

$b_k=\stackrel{\‾}{I_L\left(k\right)}-a_k{\mathrm{\μ}}_k$

其中，μ_k表示均值比值差异图I_M在ω_k邻域窗口中像素的平均值，表示对数比值差异图I_L在ω_k邻域窗口中像素的平均值，表示均值比值差异图I_M在ω_k邻域窗口中像素的方差，ε取值为10^-6。

步骤5：利用水平集模型分割导向差异图I_G。

利用基于局部相关熵的ε均值聚类水平集模型对导向差异图I_G进行分割，得到变化检测的结果，具体步骤如下：

(5a)初始化水平集函数为Φ＝Φ⁰，并根据水平集函数Φ在导向差异图I_G上做一条闭合曲线：

${\mathrm{\Φ}}^0=\mathrm{\Φ}(m,t=0)=\left\{\begin{array}{cc}-c& m\∈{\mathrm{\Ω}}_1\\ 0& m\∈{\mathrm{\Ω}}_2\\ c& m\∈{\mathrm{\Ω}}_3\end{array}\right.$

其中，m表示导向差异图I_G中像素点，c表示大于0的常数，本实例中c取常数1，Ω₁表示水平集函数Φ将导向差异图I_G分成的闭合曲线内部区域像素点，Ω₂表示水平集函数Φ将导向差异图I_G分成的闭合曲线上像素点，Ω₃表示水平集函数Φ将导向差异图I_G分成的闭合曲线外部区域像素点；

(5b)初始化权重因子其中，表示第i个像素点在第t次迭代中的权重因子，t的初始值为0，N表示导向差异图I_G的总像素点数，j表示导向差异图I_G中第j个像素点；

(5c)分别计算当前水平集函数Φ将导向差异图I_G分成的内部区域的灰度均值和外部区域的灰度均值

${\mathrm{\α}}_{j,1}^t=\frac{{\mathrm{\Σ}}_{i\∈G_j}k(j,i)H_{\mathrm{\ϵ}}\left({\mathrm{\Φ}}_j^{t-1}\right)w_i^t{x}_i}{{\mathrm{\Σ}}_{i\∈G_j}k(j,i)H_{\mathrm{\ϵ}}\left({\mathrm{\Φ}}_j^{t-1}\right)w_i^t}$

${\mathrm{\α}}_{j,2}^t=\frac{{\mathrm{\Σ}}_{i\∈G_j}k(j,i)(1-H_{\mathrm{\ϵ}}\left({\mathrm{\Φ}}_j^{t-1}\right))w_i^t{x}_i}{{\mathrm{\Σ}}_{i\∈G_j}k(j,i)\left({1-H}_{\mathrm{\ϵ}}\left({\mathrm{\Φ}}_j^{t-1}\right)\right)w_i^t}$

其中，t表示迭代次数，初始值为0，j表示导向差异图I_G中第j个像素点，G_j表示以第j个像素点为中心，大小为3×3的邻域窗口，k(j,i)为高斯核函数，H_ε(x)为Heavisirle函数，表示第j个像素点在第t-1次迭代中的水平集函数值，表示在领域窗口G_j中第i个像素点在第t次迭代中的权重因子，初始值为x_i表示在领域窗口G_j中第i个像素点的像素值。

(5d)按照下式，计算权重因子

${\left\{w_j^t\right\}}_{j=1}^N=\underset{i\∈G_j}{\mathrm{\Σ}}(H_{\mathrm{\ϵ}}\left({\mathrm{\Φ}}_j^{t-1}\right)f\left({\left|\right|x_i-{\mathrm{\α}}_{j,1}^t\left|\right|}_2\right)+(1-H_{\mathrm{\ϵ}}\left({\mathrm{\Φ}}_j^{t-1}\right))f\left({\left|\right|x_i-{\mathrm{\α}}_{j,2}^t\left|\right|}_2\right)),$

其中，N表示导向差异图I_G的总像素点数，f(x)＝exp(-x²/2σ²)表示窗口大小为σ的高斯核函数，σ＝3，||·||₂表示2范式。

(5e)通过判别规则，收缩权重因子判别规则用公式表示如下：

$w_i^t=\left\{\begin{array}{cc}{w}_i^t,& w_i^t\&GreaterEqual;\mathrm{\&tau;}\\ 0,& w_i^t<\mathrm{\&tau;}\end{array}\right.$

其中，τ表示阈值，τ取0.1。

(5f)利用如下水平集函数更新公式对当前水平集函数Φ进行更新：

${\mathrm{\&Phi;}}_j^t={\mathrm{\&Phi;}}_j^{t-1}-\mathrm{\&lambda;}\frac{\mathrm{dE}(\mathrm{\&Phi;},{\{{\mathrm{\&alpha;}}_{j,1}^t,{\mathrm{\&alpha;}}_{j,2}^t\}}_{j=1}^N)}{d{\mathrm{\&Phi;}}_j}$

$\begin{array}{c}\frac{\mathrm{dE}(\mathrm{\&Phi;},{\{{\mathrm{\&alpha;}}_{j,1}^t,{\mathrm{\&alpha;}}_{j,2}^t\}}_{j=1}^N)}{d{\mathrm{\&Phi;}}_j}=H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j\right){\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j\right)w_j^t\underset{i\&Element;G_j}{\mathrm{\&Sigma;}}k(j,i)({\left|\right|x_j-{\mathrm{\&alpha;}}_{j,1}^t\left|\right|}_2^2-{\left|\right|x_j-{\mathrm{\&alpha;}}_{j,2}^t\left|\right|}_2^2)\\ -\mathrm{\&upsi;}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j\right)\mathrm{div}\left(\frac{\&dtri;{\mathrm{\&Phi;}}_j}{\&dtri;\left|{\mathrm{\&Phi;}}_j\right|}\right)-\mathrm{\&beta;}({\&dtri;}^2{\mathrm{\&Phi;}}_j-\mathrm{div}\left(\frac{\&dtri;{\mathrm{\&Phi;}}_j}{\&dtri;\left|{\mathrm{\&Phi;}}_j\right|}\right))\end{array}$

其中，λ表示步长，本发明λ取0.1，表示梯度下降流公式，表示Dirac函数，ε取1.0，υ表示变量系数，υ取值为0.003×255²，div(·)表示散度，▽表示梯度，β表示变量系数，β取值为3.0，▽²表示拉普拉斯算子。

(5g)对步骤(5c)～(5f)进行重复迭代操作，并判断当前水平集函数Φ是否收敛即Φ^t＝Φ^t-1，如果是，则停止对水平集函数Φ的更新，得到变化检测结果图，完成SAR图像的变化检测；如果否，则继续通过循环迭代对水平集函数进行更新。

本发明的效果可以通过以下仿真做进一步说明：

1.仿真条件：

本发明的硬件测试平台是：处理器为Inter Core i3350M，主频为2.27GHz，内存2GB，软件平台为：Windows 7旗舰版32位操作系统和Matlab R2010b。本发明的输入图像分别为Bern地区不同时刻的两幅SAR图像和Ottawa地区不同时刻的两幅SAR图像，其中Bern地区不同时刻的两幅SAR图像，图像大小均为301×301，灰度级为256，格式为BMP，Ottawa地区不同时刻的两幅SAR图像，图像大小均为290×350，灰度级为256，格式为BMP。

2.仿真参数：

分别对Bern地区和Ottawa地区不同时刻的SAR图像进行仿真实验，其实验结果图利用如下参数来进行定量分析：

①计算漏检数：把参考标准图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检数FN。

②计算误检数：把参考标准图中未发生变化但实验结果图中检测为变化的像素个数，称为误检数FP。

③总误差数：计算漏检数和误检数的和，称为总误差数OE。

④正确率：

⑤Kappa系数： $\mathrm{Kappa}=\frac{\mathrm{Pcc}-\mathrm{PC}*\mathrm{pc}}{\mathrm{PC}-\mathrm{PC}*\mathrm{pc}};$

其中，Pcc表示实验结果图的变化像素和参考标准图中变化像素重合的个数占实验结果图像总像素个数的比例，PC表示参考标准图中变化像素个数占实验结果图像总像素个数的比例，pc表示实验结果图的变化像素个数占实验结果图像总像素个数的比例。

3.仿真方法：

方法1：A.Ghosh等学者在文章“Fuzzy clustering algorithms for unsupervisedchange detection in remote sensing images,Inform.Sci.,vol.181,no.4,pp.699–715,Feb.2011.”中提出的基于FCM的变化检测方法，简称DFCM方法。

方法2：S.Krinidis等学者在文章,“A robust fuzzy local information C-meansclustering algorithm,”IEEE Trans.Image Process.,vol.19,no.5,pp.1328–1337,May2010.”中提出的FLICM方法，对对数比值差异图进行分类得到变化检测结果，简称DFLICM方法。

方法3:M.Gong等学者在文章“Change detection in synthetic aperture radar imagesbased on image fusion and fuzzy clustering,IEEE Trans.Image Process.,vol.21,no.4,pp.2141–2151,Apr.2012.”中提出的基于小波融合的变化检测方法，简称DWTF方法。

方法4：本发明方法。

4.仿真内容与结果分析：

实验1：用上述四种方法对图2(a)和图2(b)所示的为Bern地区拍摄于不同时刻的SAR图像的数据集进行变化检测的仿真实验，结果如图3和图4，其中：

图2(a)为1999年4月Bern地区河流洪水爆发前的SAR图像；

图2(b)为1999年5月Bern地区河流洪水爆发后，淹没了部分城市和机场的SAR图像；

图2(c)为Bern地区数据集变化检测的标准参考图，包含1155个变化像素和89446个非变化像素；

图3(a)表示利用图2(a)和图2(b)采用均值比值法构造的差异图像；

图3(b)表示利用图2(a)和图2(b)采用对数比值法构造的差异图像；

图3(c)表示利用图2(a)和图2(b)采用本发明方法构造的导向差异图像；

图4(a)表示利用图2(a)和图2(b)采用DFCM方法进行仿真的实验结果图；

图4(b)表示利用图2(a)和图2(b)采用DFLICM方法进行仿真的实验结果图；

图4(c)表示利用图2(a)和图2(b)采用DWTF方法进行仿真的实验结果图；

图4(d)表示利用图2(a)和图2(b)采用本发明方法进行仿真的实验结果图；

对Bern地区数据集分别采用DFCM方法，DFLICM方法，DWTF方法和本发明方法进行仿真实验的对比参数如表1所示：

表1.采用不同方法对Bern地区进行仿真的结果

从图3(a)可以看出均值比值差异图展现了全面的变化信息，但因为不能有效地抑制SAR图像的乘性相干斑噪声，造成差异图中含有大量的干扰因素。

对数比值差异图像图3(b)，虽然能够大幅度的抑制SAR图像的乘性相干斑噪声，良好的避免了背景区域对变化信息的干扰，但是也同时造成了部分变化区域的细节信息丢失。

本发明构造的差异图像图3(c)根据导向滤波器的思想将均值比值差异图作为导向图，对数比值差异图作为输入图进行融合，不仅能够有效地抑制背景区域信息，减少SAR图像乘性相干斑噪声的干扰，而且能够充分的展示变化区域信息，使得对同一区域不同时刻的SAR图像之间的变化信息有了更精确地描述。

结合图4(a)、4(b)、4(c)、4(d)和表1中总误差数OE、正确率PCC和Kappa系数可以看出，本发明方法与其他变化检测方法相比得到的变化检测结果在定量分析中总误差数OE、正确率PCC和Kappa系数都有了大幅度地提升。

实验2：用上述四种方法对图5(a)和图5(b)所示的为Ottawa地区拍摄于不同时刻的SAR图像的数据集进行变化检测的仿真实验，结果如图6和图7，其中，

图5(a)为1997年5月Ottawa地区正值洪水时期的SAR图像；

图5(b)为1997年8月Ottawa地区洪水过后的SAR图像；

图5(c)为Ottawa数据集变化检测的参考标准图，包含16049个变化像素和85451个非变化像素；

图6(a)表示利用图5(a)和图5(b)采用均值比值法构造的差异图像；

图6(b)表示利用图5(a)和图5(b)采用对数比值法构造的差异图像；

图6(c)表示利用图5(a)和图5(b)采用本发明方法构造的导向差异图像；

图7(a)表示利用图5(a)和图5(b)采用DFCM方法进行仿真的结果图；

图7(b)表示利用图5(a)和图5(b)采用DFLICM方法进行仿真的结果图；

图7(c)表示利用图5(a)和图5(b)采用DWTF方法进行仿真的结果图；

图7(d)表示利用图5(a)和图5(b)采用本发明方法进行仿真的结果图；

对Ottawa地区数据集分别采用DFCM方法，DFLICM方法，DWTF方法和本发明方法进行仿真实验，得到的对比参数如表2所示：

表2.采用不同方法对Ottawa地区进行仿真的参数结果

从图6(a)、6(b)、6(c)中可以看出，本发明构造的差异图像和实验1中定量分析的一样，不仅能够有效地抑制背景区域信息，减少SAR图像乘性相干斑噪声的干扰，而且能够充分的展示变化区域信息，使得对同一区域不同时刻的SAR图像之间的变化信息有了更精确地描述。

结合图7(a)、7(b)、7(c)、7(d)和表2中总误差数OE、正确率PCC和Kappa系数可以看出，本发明在变化区域较大的情况下也可以精确地描述变化区域信息。结合实验1和实验2可以看出本发明的方法无论在变化区域较小或者变化区域较大的情况下，变化检测性能都有了显著地提升。

Claims (8)

1.一种基于导向差异图的SAR图像变化检测方法，包括如下步骤：

(1)输入同一区域t₁时刻的SAR图像X₁和t₂时刻的SAR图像X₂，并依次进行图像配准、几何校正和图像去噪预处理，获得预处理后的两幅不同时刻的SAR图像I₁，I₂；

(2)对预处理后的两幅SAR图像I₁，I₂，通过均值比值法得到均值比值差异图I_M；

(3)对预处理后的两幅SAR图像I₁，I₂，通过对数比值法构造对数比值差异图I_L：

(4)利用导向滤波器，将均值比值差异图I_M作为导向图，将对数比值差异图I_L作为输入图，并对该导向图和输入图进行融合生成导向差异图I_G；

(5)利用基于局部相关熵的K均值聚类水平集模型对导向差异图I_G进行分割，得到变化检测的结果；

(5a)初始化水平集函数为Φ＝Φ⁰，并根据水平集函数Φ在导向差异图I_G上做一条闭合曲线：

${\mathrm{\&Phi;}}^0=\mathrm{\&Phi;}(m,t=0)=\left\{\begin{array}{cc}-c& m\&Element;{\mathrm{\&Omega;}}_1\\ 0& m\&Element;{\mathrm{\&Omega;}}_2\\ c& m\&Element;{\mathrm{\&Omega;}}_3\end{array}\right.$

其中，m表示导向差异图I_G中像素点，c表示大于0的常数，c取常数1，Ω₁表示水平集函数Φ将导向差异图I_G分成的闭合曲线内部区域像素点，Ω₂表示水平集函数Φ将导向差异图I_G分成的闭合曲线上像素点，Ω₃表示水平集函数Φ将导向差异图I_G分成的闭合曲线外部区域像素点；

(5b)初始化权重因子：其中，N表示导向差异图I_G的总像素点数，j表示导向差异图I_G中第j个像素点；

(5c)分别计算当前水平集函数Φ将导向差异图I_G分成的内部区域的灰度均值和外部区域的灰度均值t表示表示迭代次数，初始值为0；

(5d)根据内部区域的灰度均值和外部区域的灰度均值通过权重因子的更新公式，更新权重因子

(5e)通过判别规则，对更新后的权重因子进行收缩；

(5f)根据水平集函数的更新公式，对当前水平集函数Φ进行更新；

(5g)对步骤(5c)～(5f)进行重复迭代操作，并判断当前水平集函数Φ是否收敛即Φ^t＝Φ^t-1，如果是，则停止对水平集函数Φ的更新，得到变化检测结果图，完成SAR图像的变化检测；如果否，则继续通过循环迭代对水平集函数进行更新。

2.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(2)中通过均值比值法得到均值比值差异图I_M，该I_M中的每个像素I_M(i,j)通过以下公式进行：

$I_M(i,j)=1-\min (\frac{{\mathrm{\&omega;}}_1(i,j)}{{\mathrm{\&omega;}}_2(i,j)},\frac{{\mathrm{\&omega;}}_2(i,j)}{{\mathrm{\&omega;}}_1(i,j)}),$

式中，ω₁(i,j)表示预处理后t₁时刻的SAR图像I₁在以坐标点为(i,j)处的像素点为中心、邻域窗口大小为5×5的局部均值，ω₂(i,j)表示预处理后t₂时刻的SAR图像I₂在以坐标点为(i,j)处的像素点为中心、邻域窗口大小为5×5的局部均值。

3.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(3)通过对数比值法构造对数比值差异图I_L，通过以下公式进行：

$I_L=\left|\log \frac{I_2+1}{I_1+1}\right|=\log (I_2+1)-\log (I_1+1)|.$

4.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(4)中融合生成导向差异图，通过以下公式进行：

$I_G=\frac{1}{\left|\mathrm{\&omega;}\right|}\underset{k\backslash i\&Element;{\mathrm{\&omega;}}_k}{\mathrm{\&Sigma;}}(a_k{I}_M+b_k)$

其中，|ω|表示ω_k窗口中含有的像素个数，k表示导向差异图I_G的第k个像素点，ω_k表示以第k个像素点为中心、大小为25×25的邻域窗口，i为ω_k邻域窗口中第i个像素点，a_k和b_k为线性系数，表示如下：

$a_k=\frac{\frac{1}{\left|\mathrm{\&omega;}\right|}{\mathrm{\&Sigma;}}_{i\&Element;{\mathrm{\&omega;}}_k}{I}_M\left(i\right)I_L\left(i\right)-{\mathrm{\&mu;}}_k\stackrel{\&OverBar;}{I_L\left(k\right)}}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}}$

$b_k=\stackrel{\&OverBar;}{I_L\left(k\right)}-a_k{\mathrm{\&mu;}}_k$

式中，μ_k表示均值比值差异图I_M在ω_k邻域窗口中像素的平均值，表示对数比值差异图I_L在ω_k邻域窗口中像素的平均值，表示均值比值差异图I_M在ω_k邻域窗口中像素的方差，ε取值为10^-6。

5.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(5c)中计算当前水平集函数Φ将导向差异图I_G分成的内部区域的灰度均值和外部区域的灰度均值通过以下公式进行：

${\mathrm{\&alpha;}}_{j,1}^t=\frac{{\mathrm{\&Sigma;}}_{i\&Element;G_j}k(j,i)H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j^{t-1}\right)w_i^t{x}_i}{{\mathrm{\&Sigma;}}_{i\&Element;G_j}k(j,i)H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j^{t-1}\right)w_i^t}$

${\mathrm{\&alpha;}}_{j,2}^t=\frac{{\mathrm{\&Sigma;}}_{i\&Element;G_j}k(j,i)(1-H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j^{t-1}\right))w_i^t{x}_i}{{\mathrm{\&Sigma;}}_{i\&Element;G_j}k(j,i)(1-H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j^{t-1}\right))w_i^t}$

其中，t表示表示迭代次数，初始值为0，j表示导向差异图I_G中第j个像素点，G_j表示以第j个像素点为中心、大小为3×3的邻域窗口，k(j,i)为高斯核函数，H_ε(x)为Heavisirle函数，表示第j个像素点在第t-1次迭代中的水平集函数值，表示在领域窗口G_j中任意像素点i在第t次迭代中的权重因子，初始值为x_i表示在领域窗口G_j中任意像素点i的像素值。

6.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(5d)的更新权重因子通过以下公式进行：

${\left\{w_j^t\right\}}_{j=1}^N=\underset{i\&Element;G_j}{\mathrm{\&Sigma;}}(H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j^{t-1}\right)f\left({\left|\right|x_i-{\mathrm{\&alpha;}}_{j,1}^t\left|\right|}_2\right)+(1-H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j^{t-1}\right))f\left({\left|\right|x_i-{\mathrm{\&alpha;}}_{j,2}^t\left|\right|}_2\right))$

其中，N表示导向差异图I_G的总像素点数，f(x)＝exp(-x²/2σ²)表示窗口大小为σ的高斯核函数，σ取值为3，||·||₂表示2范式。

7.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(5e)的判别规则，用公式表示如下：

$w_i^t=\left\{\begin{array}{cc}{w}_i^t,& w_i^t\&GreaterEqual;\mathrm{\&tau;}\\ 0,& w_i^t<\mathrm{\&tau;}\end{array}\right.$

其中，表示像素点i在第t次迭代中的权重因子，τ表示阈值，τ取0.1。

8.根据权利要求1所述的基于导向差异图的SAR图像变化检测方法，其特征在于：所述步骤(5f)中更新公式表示如下：

${\mathrm{\&Phi;}}_j^t={\mathrm{\&Phi;}}_j^{t-1}-\mathrm{\&lambda;}\frac{\mathrm{dE}(\mathrm{\&Phi;},{\{{\mathrm{\&alpha;}}_{j,1}^t,{\mathrm{\&alpha;}}_{j,2}^t\}}_{j=1}^N)}{d{\mathrm{\&Phi;}}_j},$

$\begin{array}{c}\frac{\mathrm{dE}(\mathrm{\&Phi;},{\{{\mathrm{\&alpha;}}_{j,1}^t,{\mathrm{\&alpha;}}_{j,2}^t\}}_{j=1}^N)}{d{\mathrm{\&Phi;}}_j}=H_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j\right){\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j\right)w_j^t\underset{i\&Element;G_j}{\mathrm{\&Sigma;}}k(j,i)({\left|\right|x_j-{\mathrm{\&alpha;}}_{j,1}^t\left|\right|}_2^2-{\left|\right|x_j-{\mathrm{\&alpha;}}_{j,2}^t\left|\right|}_2^2)\\ -\mathrm{\&upsi;}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left({\mathrm{\&Phi;}}_j\right)\mathrm{div}\left(\frac{{\&dtri;\mathrm{\&Phi;}}_j}{\&dtri;\left|{\mathrm{\&Phi;}}_j\right|}\right)-\mathrm{\&beta;}({\&dtri;}^2{\mathrm{\&Phi;}}_j-\mathrm{div}(\frac{{\&dtri;\mathrm{\&Phi;}}_j}{\&dtri;\left|{\mathrm{\&Phi;}}_j\right|}\left)\right)\end{array}$

其中，λ表示步长，λ＝0.1，表示梯度下降流公式，表示Dirac函数，ε＝1.0，υ表示变量系数，υ＝0.003×255²，div(·)表示散度，▽表示梯度，β表示变量系数，β＝3.0，▽²表示拉普拉斯算子。