CN103839257A - 一种广义高斯k&i的sar图像变化检测方法 - Google Patents

Abstract

本发明提出了一种广义高斯K&I的SAR图像变化检测方法，其实现过程为：（1）首先对两幅配准后的同一地区不同时相的SAR图像采用改进的邻域比操作算子做出差异图，在对邻域处理时，与中心像素相邻的像素不再是被同等对待，而是根据它们与中心像素的相似性被赋予不同的权重；（2）通过采用一种修正的基于广义高斯模型的K&I阈值分割方法对所得的差异图进行二分类从而得到变化部分和未变化部分。本发明通过对原始K&I方法的指标函数加入一个局部能量修正项使其对邻域比差异图能找到一个合理的阈值。本发明在一定程度上改善了噪声的敏感问题，同时也显著地提高了检测效果及检测的正确率。

Description

一种广义高斯K&I的SAR图像变化检测方法

技术领域

本发明属于遥感图像处理领域，涉及遥感图像的变化检测，具体提供一种广义高斯K&I的SAR图像变化检测方法，可用于对SAR图像变化的检测。 

背景技术

合成孔径雷达（synthetic aperture radar,SAR）具有分辨率高，全天时、全天候工作的特点，与可见光、红外传感器比较具有独特的优势和无法替代的作用，被广泛应用于工农业生产、科研和军事等领域。遥感图像变化检测是通过对不同时期的两幅遥感图像的比较分析,根据图像之间的差异来得到所需的地物变化信息。随着遥感技术的发展,变化检测技术已经成为遥感图像应用的一个重要分支。目前,变化检测技术已广泛地应用于多个领域,如森林覆盖变化、土地利用覆盖、城市环境变化等,尤其是在自然灾害评估分析的应用方面。 

在变化检测研究中,国内外学者们做了大量的研究,目前SAR 图像变化检测方法可以分为:(1)基于简单代数运算的变化检测，经典的变化检测方法包括图像差值法、图像比值法、对数比值法；(2)基于图像变换的变换检测，经典的变化检测方法包括主成份分析、变化向量分析法、相关分析法图像变换法；(3)基于图像分类的变化检测方法。当前较为流行的方法是先对两幅配准图像做差异图，再对所得的差异图进行分类比较。这类方法的研究空间比较大，思路简单明确，检测精度较高。 

图像差值法：其主要过程是将两幅不同时间相同地域的遥感图像对应像素的灰度值相减得到差异影像图。理论上，在得到的差值图像上，差值为0或接 近0的认为是不变区域，不为0的认为是变化区域。这种方法的优点在于理论相对简单，容易理解和掌握。但是缺点也比较明显，由于该方法主要通过分析地物光谱值改变的灰度差值图像来实现变化检测，但在某些情况下，仅仅利用光谱特征的差值图像难以实现地物的变化检测，易产生“伪变化”信息。 

图像比值法：其主要过程是计算两幅不同时间相同地域的遥感图像对应像素灰度值的比值以得到差异影像图，其中像素比值为1或近似为1的认为是未发生变化的区域，反之，远大于或远小于1的认为是发生变化的区域。该方法可以克服乘性噪声的干扰，但该方法是假设比值图像呈正态分布的，在很多实际问题中该假设并不总是成立的，缺少适应性。 

主分量分析法：该方法又称K-L变换法，是一种经典的数学变换方法，主要过程是把原来多个波段中的有用信息集中到互不相关的新成分图像中，达到用于压缩和信息集中的目的。在进行主成分变换时，利用协方差矩阵和相关矩阵得到的主分量是不同的，由相关矩阵推导的主分量变换对于多时相分析尤其有用。 

变化向量分析法：该方法是一种多变量分析方法，描述了从第一时间到第二时间的光谱变化的强度及方向。如果变化向量的幅值超过给定的门限，则判定该像素发生变化，变化向量的方向包含变化类信息。该方法可用在多通道极化SAR图像变化检测，或者用于多特征，如空间结构特征，纹理特征分析等。 

分类后比较法：该方法是一种比较直观的方法，其优点是可以回避基于差异影像分析方法所要求的影像系列时相一致性条件,以及影像间辐射校正、匹配等问题。同时,可确定变化信息的空间位置,及提供变化信息的类型,但此方法很难获得不同时相图像的不变信息具有相同类别的分类结果,且会夸大变化 程度。同时,由于分类累积误差问题降低了变化检测精度。 

发明内容

本发明的目的提供一种利用改进邻域比及修正的一种广义高斯K&I的SAR图像变化检测方法，以便降低噪声的影响，对于服从非对称的高斯分布的邻域比图像能够找到合适的阈值，从而提高了变化检测精度。 

本发明的技术方案是，广义高斯K&I的SAR图像变化检测方法，其特征是：包括如下步骤： 

步骤101：开始基于改进邻域比及广义高斯K&I的SAR图像变化检测方法； 

步骤102：导入两幅配准后的同一地区不同时刻的SAR图像； 

步骤103：对邻域比值算法做改进，按照如下步骤构造两幅不同时间同一地区的图像I₁,I₂的差异影像图DI； 

步骤104：使用修正的广义高斯K&I算法对上述所得的差异图DI寻找最优的分割阈值T^*，由于原始的K&I阈值方法只能对对称的高斯图像找到合适的分割阈值，而通过邻域比算子所得到的图像符合非对称的高斯模型，因此对原始的广义高斯K&I方法的指标函数J_{GG_KI}(T)加入一个局部能量修正项使其对邻域比差异图能找到一个合理的阈值； 

步骤105：利用步骤104中求得的最优阈值T^*对差异图DI进行阈值分割，得到最终的变化检测结果图CNI； 

步骤106：结束基于改进邻域比及广义高斯K&I的SAR图像变化检测方法。 

所述的步骤103，包括如下步骤： 

步骤201：开始改进的邻域比算子做差异图； 

步骤202：分别取得两时像SAR图像I₁,I₂在同一位置x上的像素邻域集合 Ω₁(x)和Ω₂(x)，其大小为N×N，N∈{3,5,7,9}； 

步骤203：差异图DI在位置x上的像素灰度值DI(x)： 

$\mathrm{DI}\left(x\right)=\∂\×\frac{\min \{I_1\left(x\right),I_2\left(x\right)\}}{\max \{I_1\left(x\right),I_2\left(x\right)\}}$

其中，I₁(i)和I₂(i)分别表示图像I₁,I₂在位置x上邻域集合Ω₁(x)和Ω₂(x)的第i个元素，DI(x)值越小，说明图像I₁,I₂在位置x上的差异性越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；求解系数

其表示比值图像上像素x的邻域异质性；

表示极小值图像上位置x邻域集合的第i个元素的权重，

表示极大值图像上位置x邻域集合的第i个元素的权重,对系数

的求解采用公式： 

$\∂=\frac{\mathrm{\σ}\left(x\right)}{\mathrm{\μ}\left(x\right)}---\left(2\right)$

σ(x)和μ(x)分别表示比值图像上像素x邻域的方差和均值。

值越大，表示区域异质性越强，差异图DI(x)的求解公式中第一部分所占比重大；反之，则第二部分所占比重大； 

步骤204：求解邻域比部分中邻域内第i个元素的权重α_i,x，它取决于邻域内像素i和中心像素x的相似性，该相似性采用加权欧式距离测量v(Ω_k)表示中心像素为k的邻域矢量。α_i,x计算公式如下： 

${\mathrm{\α}}_{i,x}=\frac{1}{Z\left(i\right)}{e}^{\frac{-{\left|\right|v\left({\mathrm{\Ω}}_i\right)-v\left({\mathrm{\Ω}}_x\right)\left|\right|}_2^2}{h^2}}---\left(3\right)$

其中Z(i)为归一化常数， 

$Z\left(i\right)=\underset{x}{\mathrm{\Σ}}{e}^{\frac{-{\left|\right|v\left({\mathrm{\Ω}}_i\right)-v\left({\mathrm{\Ω}}_x\right)\left|\right|}_2^2}{h^2}}---\left(4\right)$

参数h的评估使用如下公式： 

$h=\sqrt{\frac{{2l}^2}{\ln (1/\mathrm{\γ})}}{\mathrm{\σ}}_n---\left(5\right)$

其中σ_n是噪声的标准差，可以直接从图像中计算出来，l²是领域大小，0＜γ≤1是自由参数，通常取0.6-0.9； 

步骤205：比较两个邻域集合Ω₁(x)和Ω₂(x)的相似性，将步骤203和步骤204中公式(2)和(3)计算的结果代入公式(1)，对图像I₁,I₂的每一个像素点计算得到差异图DI； 

步骤206：结束改进的邻域比算子做差异图。 

所述的步骤104，包括如下步骤： 

步骤301：开始修正的广义高斯K&I方法求图像阈值； 

步骤302：求得差异图DI的灰度直方图h(X_l)，l∈{0,1,......,L-1}； 

步骤303：基于上述所得的灰度直方图h(X_l)，对每一个阈值T(T∈{0,1,......,L-1})，灰度值小于T的属于变化类，灰度值大于T的属于未变化类，分别求解变化类和未变化类的概率P_u(T)和P_c(T)，均值m_u(T)和m_c(T)以及方差

和

$\left\{\begin{array}{c}{P}_u\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right),m_u\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right)\\ {\mathrm{\σ}}_u^2\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_{l=0}}{\overset{T}{\mathrm{\Σ}}}{[X_l-m_u\left(T\right)]}^{2}h\left(X_l\right)\\ P_c\left(T\right)=1-P_u\left(T\right),m_c\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right)\\ {\mathrm{\σ}}_c^2\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{[X_l-m_c\left(T\right)]}^{2}h\left(X_l\right)\end{array}\right.---\left(6\right)$

步骤304：使用查表法求解广义高斯建模下两类条件概率表达式  $P\left(X_l\right|{\mathrm{\ω}}_i)=a_i{e}^{-{[b_i{|X}_l-m_i|]}^{{\mathrm{\β}}_i}}$ 的形状参数 $a_i=\frac{b_i{\mathrm{\β}}_i}{2\mathrm{\Γ}\left(\frac{1}{{\mathrm{\β}}_i}\right)}$ 和 $b_i=\frac{1}{{\mathrm{\σ}}_i}\sqrt{\frac{\mathrm{\Γ}\left(\frac{3}{{\mathrm{\β}}_i}\right)}{\mathrm{\Γ}\left(\frac{1}{{\mathrm{\β}}_i}\right)}}(i=u,c),$ 其中β_i为分布的形状参数； 

步骤305：根据步骤303和步骤304的计算结果求解原始的广义高斯K&I的指标函数J_{GG_KI}(T)： 

$J_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right){\left[b_u\left(T\right)\right|X_l-m_u\left(T\right)\left|\right]}^{{\mathrm{\β}}_u\left(T\right)}$

$+\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}h\left(X_l\right){\left[b_u\left(T\right)\right|X_l-m_c\left(T\right)\left|\right]}^{{\mathrm{\β}}_u\left(T\right)}$

$+H(\mathrm{\ω},T)-[P_u\left(T\right)\ln a_u\left(T\right)+P_c\left(T\right)\ln a_c\left(T\right)]---\left(7\right)$

步骤306：计算利用了空间信息的局部能量修正项

其中S是图像的所有像素集合，定义

A_i(i∈S)代表变化类或未变化类的标签(像素i属于为变化类，A_i＝0；否则，A_i＝1。如果A_ξ＝A_η，[A_ξ≠A_η]为0，如果A_ξ≠A_η,[A_ξ≠A_η]为1)。系数B(ξ,η）=exp(-θ||Y_ξ-Y_η||²)是对像素ξ和η之 间不连续性的惩罚,ξ和η越相似，B(ξ,η）值越大，反之，B(ξ,η）值接近于0，这里θ≥0是平滑参数，Y_i(i∈S)表示像素i的灰度值； 

步骤307：对每一个阈值T(T∈{0,1,......,L-1})，计算修正后的指标函数  $J\left(T\right)={\mathrm{\λJ}}_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)+\underset{\mathrm{\η}\∈S}{\mathrm{\Σ}}{B}_{\mathrm{\η}}\left(T\right);$

步骤308：根据 $T^{*}=\arg \underset{T=\mathrm{0,1},......,L-1}{\min J\left(T\right)}$ 最小化指标函数 $J\left(T\right)=\mathrm{\λ}{J}_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)+\underset{\mathrm{\η}\∈S}{\mathrm{\Σ}}{B}_{\mathrm{\η}}\left(T\right),$ 从而求得最优分割阈值T^*； 

步骤309：结束修正的广义高斯KI算法求分割阈值。 

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

1、本发明是在邻域比的基础上进行改进构造差异影像图，在对邻域处理时邻域内的每一个像素都被赋予了一个合理的权重，更加有效的利用了空间信息，进一步降低了变化检测过程中斑点噪声的影响； 

2、本发明针对改进的邻域比方法构造的差异影像图，提出了一种修正的广义高斯K&I方法。由于原始的广义高斯K&I方法不适合这种符合非对称的高斯模型的差异图，通过局部能量项的修正，使其能够很好的找出优良的分割阈值； 

3、仿真结果表明，本发明采用的改进邻域比及修正的广义高斯K&I方法与对数比值LR算子，原始的邻域比NR算子构造差异影像图方法及广义K&I（GKI），广义高斯K&I（GG_KI）方法相比，正确检测率高，Kappa系数高。 

附图说明

图1是基于改进邻域比及广义高斯K&I的SAR图像变化检测方法主流程图； 

图2是改进邻域比算子的流程图； 

图3是修正广义高斯K&I算法的流程图； 

图4是第一组实验仿真图，图4（a）和图4（b）的拍摄时间分别为1999.04 和1999.05，大小均为301×301，图4（c）为参考图； 

图5是对应图4的不同方法构造的差异影像图，其中图5（a），图5（b）和图5（c）分别表示由LR算子，NR算子和本发明方法WNR构造的差异影像图； 

图6是对应图5的不同方法构造的差异影像图的ROC曲线的比较； 

图7是对应图5的GKI阈值分割变化检测结果图，其中图7（a），图7（b）图7（c）分别表示由LR算子，NR算子和本发明方法WNR构造的差异影像图在经过GKI阈值分割结果图； 

图8是第二组实验仿真图，图8（a）和图8（b）的拍摄时间分别为1997.05和1997.08，大小均为290×350，图8（c）为参考图； 

图9是对应图8的基于的变化检测结果图，其中图9（a），图9（b）和图9（c）分别表示由GKI方法，GG_KI方法和本发明方法即修正的广义高斯K&I方法处理后得到的阈值分割变化检测结果图； 

图10是基于改进邻域比差异图所得到的变化检测结果与参考图的直方图比较，图10（a），图10（b）分别是用GG_KI方法和本发明方法做差异图分析。 

具体实施方式

一种广义高斯K&I的SAR图像变化检测方法，关键是要找到合适的图像分割阈值。由于原始的K&I阈值方法只能对对称的高斯图像分割找到合适的阈值，而通过邻域比算子所得到的图像符合非对称的高斯模型，因此对原始K&I方法的指标函数加入一个局部能量修正项使其对邻域比差异图能找到一个合理的阈 值。本发明应用于SAR图像的变化检测，一定程度上改善了噪声的敏感问题，同时也显著地提高了检测效果及检测的正确率。 

一种广义高斯K&I的SAR图像变化检测方法的特征是：首先对两幅配准后的同一地区不同时相的SAR图像采用改进的邻域比操作算子做出差异图，在对邻域处理时，与中心像素相邻的像素不再是被同等对待，而是根据它们与中心像素的相似性被赋予不同的权重；（2）通过采用一种修正的基于广义高斯模型的K&I阈值分割方法对所得的差异图进行二分类从而得到变化部分和未变化部分。 

如图1所示。 

主流程图步骤特征是： 

步骤101：开始基于改进邻域比及广义高斯K&I的SAR图像变化检测方法； 

步骤102：导入两幅配准后的同一地区不同时刻的SAR图像； 

步骤103：对邻域比值算法做改进，按照如下步骤构造两幅不同时间同一地区的图像I₁,I₂的差异影像图DI； 

步骤104：使用修正的广义高斯K&I算法对上述所得的差异图DI寻找最优的分割阈值T^*，由于原始的K&I阈值方法只能对对称的高斯图像找到合适的分割阈值，而通过邻域比算子所得到的图像符合非对称的高斯模型，因此对原始的广义高斯K&I方法的指标函数J_{GG_KI}(T)加入一个局部能量修正项使其对邻域比差异图能找到一个合理的阈值； 

步骤105：利用步骤104中求得的最优阈值T^*对差异图DI进行阈值分割，得到最终的变化检测结果图CNI； 

步骤106：结束基于改进邻域比及广义高斯K&I的SAR图像变化检测方法。 

如图2所示。 

所述的步骤103，包括如下步骤： 

步骤201：开始改进的邻域比算子做差异图； 

步骤202：分别取得两时像SAR图像I₁,I₂在同一位置x上的像素邻域集合Ω₁(x)和Ω₂(x)，其大小为N×N，N∈{3,5,7,9}； 

步骤203：差异图DI在位置x上的像素灰度值DI(x)： 

$\mathrm{DI}\left(x\right)=\∂\×\frac{\min \{I_1\left(x\right),I_2\left(x\right)\}}{\max \{I_1\left(x\right),I_2\left(x\right)\}}$

其中，I₁(i)和I₂(i)分别表示图像I₁,I₂在位置x上邻域集合Ω₁(x)和Ω₂(x)的第i个元素，DI(x)值越小，说明图像I₁,I₂在位置x上的差异性越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；求解系数

其表示比值图像上像素x的邻域异质性；

表示极小值图像上位置x邻域集合的第i个元素的权重，

表示极大值图像上位置x邻域集合的第i个元素的权重,对系数

的求解采用公式： 

$\∂=\frac{\mathrm{\σ}\left(x\right)}{\mathrm{\μ}\left(x\right)}---\left(2\right)$

σ(x)和μ(x)分别表示比值图像上像素x邻域的方差和均值。

值越大，表示区域异质性越强，差异图DI(x)的求解公式中第一部分所占比重大；反之，则第二部分所占比重大； 

步骤204：求解邻域比部分中邻域内第i个元素的权重α_i,x，它取决于邻域内像素i和中心像素x的相似性，该相似性采用加权欧式距离测量 v(Ω_k)表示中心像素为k的邻域矢量。α_i,x计算公式如下： 

${\mathrm{\α}}_{i,x}=\frac{1}{Z\left(i\right)}{e}^{\frac{-{\left|\right|v\left({\mathrm{\Ω}}_i\right)-v\left({\mathrm{\Ω}}_x\right)\left|\right|}_2^2}{h^2}}---\left(3\right)$

其中Z(i)为归一化常数， 

$Z\left(i\right)=\underset{x}{\mathrm{\Σ}}{e}^{\frac{-{\left|\right|v\left({\mathrm{\Ω}}_i\right)-v\left({\mathrm{\Ω}}_x\right)\left|\right|}_2^2}{h^2}}---\left(4\right)$

参数h的评估使用如下公式： 

$h=\sqrt{\frac{{2l}^2}{\ln (1/\mathrm{\γ})}}{\mathrm{\σ}}_n---\left(5\right)$

其中σ_n是噪声的标准差，可以直接从图像中计算出来，l²是领域大小，0＜γ≤1是自由参数，通常取0.6-0.9； 

步骤205：比较两个邻域集合Ω₁(x)和Ω₂(x)的相似性，将步骤203和步骤204中公式(2)和(3)计算的结果代入公式(1)，对图像I₁,I₂的每一个像素点计算得到差异图DI； 

步骤206：结束改进的邻域比算子做差异图。 

如图3所示。 

所述的步骤104，包括如下步骤： 

步骤301：开始修正的广义高斯K&I方法求图像阈值； 

步骤302：求得差异图DI的灰度直方图h(X_l)，l∈{0,1,......,L-1}； 

步骤303：基于上述所得的灰度直方图h(X_l)，对每一个阈值T(T∈{0,1,......,L-1})，灰度值小于T的属于变化类，灰度值大于T的属于未变化类，分别求解变化类和未变化类的概率P_u(T)和P_c(T)，均值m_u(T)和m_c(T)以及方差和

$\left\{\begin{array}{c}{P}_u\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right),m_u\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right)\\ {\mathrm{\σ}}_u^2\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_{l=0}}{\overset{T}{\mathrm{\Σ}}}{[X_l-m_u\left(T\right)]}^{2}h\left(X_l\right)\\ P_c\left(T\right)=1-P_u\left(T\right),m_c\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right)\\ {\mathrm{\σ}}_c^2\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{[X_l-m_c\left(T\right)]}^{2}h\left(X_l\right)\end{array}\right.---\left(6\right)$

步骤304：使用查表法求解广义高斯建模下两类条件概率表达式  $P\left(X_l\right|{\mathrm{\ω}}_i)=a_i{e}^{-{[b_i{|X}_l-m_i|]}^{{\mathrm{\β}}_i}}$ 的形状参数 $a_i=\frac{b_i{\mathrm{\β}}_i}{2\mathrm{\Γ}\left(\frac{1}{{\mathrm{\β}}_i}\right)}$ 和 $b_i=\frac{1}{{\mathrm{\σ}}_i}\sqrt{\frac{\mathrm{\Γ}\left(\frac{3}{{\mathrm{\β}}_i}\right)}{\mathrm{\Γ}\left(\frac{1}{{\mathrm{\β}}_i}\right)}}(i=u,c),$ 其中β_i为分布的形状参数； 

步骤305：根据步骤303和步骤304的计算结果求解原始的广义高斯K&I的指标函数J_{GG_KI}(T)： 

$J_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right){\left[b_u\left(T\right)\right|X_l-m_u\left(T\right)\left|\right]}^{{\mathrm{\β}}_u\left(T\right)}$

$+\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}h\left(X_l\right){\left[b_u\left(T\right)\right|X_l-m_c\left(T\right)\left|\right]}^{{\mathrm{\β}}_u\left(T\right)}$

$+H(\mathrm{\ω},T)-[P_u\left(T\right)\ln a_u\left(T\right)+P_c\left(T\right)\ln a_c\left(T\right)]---\left(7\right)$

步骤306：计算利用了空间信息的局部能量修正项

其中S是图像的所有像素集合，定义

A_i(i∈S)代表变化类或未变化类的标签(像素i属于为变化类，A_i＝0；否则，A_i＝1。如果A_ξ＝A_η，[A_ξ≠A_η]为0，如果A_ξ≠A_η,[A_ξ≠A_η]为1)。系数B(ξ,η）=exp(-θ||Y_ξ-Y_η||²)是对像素ξ和η之 间不连续性的惩罚,ξ和η越相似，B(ξ,η）值越大，反之，B(ξ,η）值接近于0，这里θ≥0是平滑参数，Y_i(i∈S)表示像素i的灰度值； 

步骤307：对每一个阈值T(T∈{0,1,......,L-1})，计算修正后的指标函数  $J\left(T\right)={\mathrm{\λJ}}_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)+\underset{\mathrm{\η}\∈S}{\mathrm{\Σ}}{B}_{\mathrm{\η}}\left(T\right);$

步骤308：根据 $T^{*}=\arg \underset{T=\mathrm{0,1},......,L-1}{\min J\left(T\right)}$ 最小化指标函数 $J\left(T\right)={\mathrm{\λJ}}_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)+\underset{\mathrm{\η}\∈S}{\mathrm{\Σ}}{B}_{\mathrm{\η}}\left(T\right),$ 从而求得最优分割阈值T^*； 

步骤309：结束修正的广义高斯KI算法求分割阈值。 

本实施实例没有详细叙述的部分属行业公知的常用手段，这里不一一叙述。 

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

1、仿真参数 

对于具有参考图的实验仿真图，可进行定量的变化检测结果分析： 

①使用受试者工作特征曲线ROC来衡量差异图的性能，真正类率TPR和负正类率FPR以及曲线面积AUC来定量衡量。 

②计算漏检个数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检个数FN； 

③计算错检个数：统计实验结果图中未发生变化区域的像素个数，与参考图中未变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为错检个数FP； 

④正确分类的概率PCC：PCC=(TP+TN)/(TP+FP+TN+FN)； 

⑤衡量检测结果图与参考图一致性的Kappa系数：Kappa=(PCC-PRE)/(1-PRE),其中 

$\mathrm{PRE}=\frac{(\mathrm{TP}+\mathrm{FP})\·\mathrm{Nc}+(\mathrm{FN}+\mathrm{TN})\·\mathrm{Nu}}{N^2}.$

2、仿真内容 

①将本发明方法改进邻域比WNR算子构造的差异影像图，与对数比值LR算子，原始邻域比值NR算子，构造的差异影像图的ROC曲线指标对比，以及使用GKI变化检测方法的结果对比和定量分析通过第一组实验仿真图完成； 

②对于分析差异图，将本发明提出的修正的广义高斯K&I方法与广义K&I（GKI），广义高斯K&I（GG_KI）的变化检测效果对比和定量分析，通过第二组实验仿真图完成。 

3、仿真实验结果及分析 

①反映Bern城市的SAR图像如图4所示，图4（a）和图4（b）的拍摄时间分别为1999.04和1999.05，大小均为301×301，图4（c）为参考图。 

用不同算子构造的差异影像图如图5所示，其中图5（a），图5（b）和图5（c）分别表示由LR算子，NR算子和本发明方法WNR构造的差异影像图。图6给出了三种方法构造所得差异影像图的ROC曲线的比较，从图中可以看出，WNR方法得到的差异图优于LR方法和NR方法的结果。用不同算子构造的差异影像图在经过GKI方法阈值分割结果如图7所示，其中图7（a），图7（b）和图7（c）分别表示由LR算子，NR算子和本发明方法构造的差异影像图在经过GKI方法阈值分割结果图，从图7可以看出，本发明方法构造的差异影像图的GKI方法结果图噪声点较少，细节清晰，变化检测效果比较好。用不同算子构造差异影像图的ROC曲线指标和经过GKI方法的变化检测结果分析分别如表1和表2所示。 

表1 不同算子构造差异影像图的ROC曲线指标的比较 

从表1中可以看出，本发明方法所得差异图的ROC曲线面积较大，负正类率FPR值较小，说明该差异图与参考图更接近，比起LR差异图和NR差异图更有优势。 

表2 不同算子构造差异影像图的变化检测结果比较 

从表2中可以看出，本发明方法的变化检测正确率以及Kappa值，比其他对比实验方法都要高，而且在误检数与漏检数上比其他两种方法都低。同时与参考图像图7(c)对比可看出，本发明对杂点的去除效果比其他两个方法都要好。 

②反映Ottawa地区水灾的SAR图像如图8所示，图8（a）和图8（b）的拍摄时间分别为1997.05和1997.08，大小均为290×350，图9（c）为参考图，图8（d）为使用WNR方法得到的差异图。 

采用WNR方法所得差异影像图经过不同的差异图分析方法所得的变化检测结果图如图9所示，其中图9（a），图9（b）和图9（c）分别表示采用GKI， GG_KI以及本发明方法所得的变化检测结果图。通过与参考图像图8（c）对比，可以看出本发明方法无论是对噪点的鲁棒性还是对细节的保持都比其他两种方法要好。图10给出了GG_KI方法和本发明方法对WNR差异图的处理结果，图10（a）是GG_KI方法的结果，图10（b）是本发明方法的结果，显然本发明方法能够找到更为合理的阈值。用不同阈值方法得到的变化检测结果分析如表3所示。 

表3 不同算法得到的变化检测结果分析 

从表3可以看出本发明方法找到的分割阈值更为合理，变化检测结果的正确率和Kappa值都高于另外两种方法，并且漏检数FN和错检数FP相对较均衡，得到的变化检测结果图比较稳定。 

Claims (3)

1.一种广义高斯K&I的SAR图像变化检测方法，其特征是：包括如下步骤：

步骤101：开始基于改进邻域比及广义高斯K&I的SAR图像变化检测方法；

步骤102：导入两幅配准后的同一地区不同时刻的SAR图像；

步骤103：对邻域比值算法做改进，按照如下步骤构造两幅不同时间同一地区的图像I₁,I₂的差异影像图DI；

步骤104：使用修正的广义高斯K&I算法对上述所得的差异图DI寻找最优的分割阈值T^*，由于原始的K&I阈值方法只能对对称的高斯图像找到合适的分割阈值，而通过邻域比算子所得到的图像符合非对称的高斯模型，因此对原始的广义高斯K&I方法的指标函数J_{GG_KI}(T)加入一个局部能量修正项使其对邻域比差异图能找到一个合理的阈值；

步骤105：利用步骤104中求得的最优阈值T^*对差异图DI进行阈值分割，得到最终的变化检测结果图CNI；

步骤106：结束基于改进邻域比及广义高斯K&I的SAR图像变化检测方法。

2.根据权利要求1所述的一种广义高斯K&I的SAR图像变化检测方法，其特征是：所述的步骤103，包括如下步骤：

步骤201：开始改进的邻域比算子做差异图；

步骤202：分别取得两时像SAR图像I₁,I₂在同一位置x上的像素邻域集合Ω₁(x)和Ω₂(x)，其大小为N×N，N∈{3,5,7,9}；

步骤203：差异图DI在位置x上的像素灰度值DI(x)：

$\mathrm{DI}\left(x\right)=\∂\×\frac{\min \{I_1\left(x\right),I_2\left(x\right)\}}{\max \{I_1\left(x\right),I_2\left(x\right)\}}$

其中，I₁(i)和I₂(i)分别表示图像I₁,I₂在位置x上邻域集合Ω₁(x)和Ω₂(x)的第i个元素，DI(x)值越小，说明图像I₁,I₂在位置x上的差异性越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；求解系数

其表示比值图像上像素x的邻域异质性；

表示极小值图像上位置x邻域集合的第i个元素的权重，表示极大值图像上位置x邻域集合的第i个元素的权重,对系数

的求解采用公式：

$\∂=\frac{\mathrm{\σ}\left(x\right)}{\mathrm{\μ}\left(x\right)}---\left(2\right)$

σ(x)和μ(x)分别表示比值图像上像素x邻域的方差和均值。

值越大，表示区域异质性越强，差异图DI(x)的求解公式中第一部分所占比重大；反之，则第二部分所占比重大；

步骤204：求解邻域比部分中邻域内第i个元素的权重α_i,x，它取决于邻域内像素i和中心像素x的相似性，该相似性采用加权欧式距离测量

v(Ω_k)表示中心像素为k的邻域矢量。α_i,x计算公式如下：

${\mathrm{\α}}_{i,x}=\frac{1}{Z\left(i\right)}{e}^{\frac{-{\left|\right|v\left({\mathrm{\Ω}}_i\right)-v\left({\mathrm{\Ω}}_x\right)\left|\right|}_2^2}{h^2}}---\left(3\right)$

其中Z(i)为归一化常数，

$Z\left(i\right)=\underset{x}{\mathrm{\Σ}}{e}^{\frac{-{\left|\right|v\left({\mathrm{\Ω}}_i\right)-v\left({\mathrm{\Ω}}_x\right)\left|\right|}_2^2}{h^2}}---\left(4\right)$

参数h的评估使用如下公式：

$h=\sqrt{\frac{{2l}^2}{\ln (1/\mathrm{\γ})}}{\mathrm{\σ}}_n---\left(5\right)$

其中σ_n是噪声的标准差，可以直接从图像中计算出来，l²是领域大小，0＜γ≤1是自由参数，通常取0.6-0.9；

步骤205：比较两个邻域集合Ω₁(x)和Ω₂(x)的相似性，将步骤203和步骤204中公式(2)和(3)计算的结果代入公式(1)，对图像I₁,I₂的每一个像素点计算得到差异图DI；

步骤206：结束改进的邻域比算子做差异图。

3.根据权利要求1所述的一种广义高斯K&I的SAR图像变化检测方法，其特征是：所述的步骤104，包括如下步骤：

步骤301：开始修正的广义高斯K&I方法求图像阈值；

步骤302：求得差异图DI的灰度直方图h(X_l)，l∈{0,1,......,L-1}；

步骤303：基于上述所得的灰度直方图h(X_l)，对每一个阈值T(T∈{0,1,......,L-1})，灰度值小于T的属于变化类，灰度值大于T的属于未变化类，分别求解变化类和未变化类的概率P_u(T)和P_c(T)，均值m_u(T)和m_c(T)以及方差

和

$\left\{\begin{array}{c}{P}_u\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right),m_u\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right)\\ {\mathrm{\σ}}_u^2\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_{l=0}}{\overset{T}{\mathrm{\Σ}}}{[X_l-m_u\left(T\right)]}^{2}h\left(X_l\right)\\ P_c\left(T\right)=1-P_u\left(T\right),m_c\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right)\\ {\mathrm{\σ}}_c^2\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{[X_l-m_c\left(T\right)]}^{2}h\left(X_l\right)\end{array}\right.---\left(6\right)$

步骤304：使用查表法求解广义高斯建模下两类条件概率表达式 $P\left(X_l\right|{\mathrm{\ω}}_i)=a_i{e}^{-{[b_i{|X}_l-m_i|]}^{{\mathrm{\β}}_i}}$ 的形状参数 $a_i=\frac{b_i{\mathrm{\β}}_i}{2\mathrm{\Γ}\left(\frac{1}{{\mathrm{\β}}_i}\right)}$ 和 $b_i=\frac{1}{{\mathrm{\σ}}_i}\sqrt{\frac{\mathrm{\Γ}\left(\frac{3}{{\mathrm{\β}}_i}\right)}{\mathrm{\Γ}\left(\frac{1}{{\mathrm{\β}}_i}\right)}}(i=u,c),$ 其中β_i为分布的形状参数；

步骤305：根据步骤303和步骤304的计算结果求解原始的广义高斯K&I的指标函数J_{GG_KI}(T)：

$J_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right){\left[b_u\left(T\right)\right|X_l-m_u\left(T\right)\left|\right]}^{{\mathrm{\β}}_u\left(T\right)}$

$+\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}h\left(X_l\right){\left[b_u\left(T\right)\right|X_l-m_c\left(T\right)\left|\right]}^{{\mathrm{\β}}_u\left(T\right)}$

$+H(\mathrm{\ω},T)-[P_u\left(T\right)\ln a_u\left(T\right)+P_c\left(T\right)\ln a_c\left(T\right)]---\left(7\right)$

步骤306：计算利用了空间信息的局部能量修正项其中S是图像的所有像素集合，定义

A_i(i∈S)代表变化类或未变化类的标签(像素i属于为变化类，A_i＝0；否则，A_i＝1。如果A_ξ＝A_η，[A_ξ≠A_η]为0，如果A_ξ≠A_η,[A_ξ≠A_η]为1)。系数B(ξ,η）=exp(-θ||Y_ξ-Y_η||²)是对像素ξ和η之间不连续性的惩罚,ξ和η越相似，B(ξ,η）值越大，反之，B(ξ,η）值接近于0，这里θ≥0是平滑参数，Y_i(i∈S)表示像素i的灰度值；

步骤307：对每一个阈值T(T∈{0,1,......,L-1})，计算修正后的指标函数 $J\left(T\right)={\mathrm{\λJ}}_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)+\underset{\mathrm{\η}\∈S}{\mathrm{\Σ}}{B}_{\mathrm{\η}}\left(T\right);$

步骤308：根据 $T^{*}=\arg \underset{T=\mathrm{0,1},......,L-1}{\min J\left(T\right)}$ 最小化指标函数 $J\left(T\right)=\mathrm{\λ}{J}_{\mathrm{GG}\_\mathrm{KI}}\left(T\right)+\underset{\mathrm{\η}\∈S}{\mathrm{\Σ}}{B}_{\mathrm{\η}}\left(T\right),$ 从而求得最优分割阈值T^*；

步骤309：结束修正的广义高斯KI算法求分割阈值。

Images