基于广义Gamma模型的KL变化信息提取及非监督KI阈值分割的变化

基金项目：国家863项目（No.2009AA12z118）、国家自然科学基金（No.40601058,40701108,40871191）作者简介：高丛珊（1985—），女，硕士研究生，现就读于中国科学院对地观测与数字地球科学中心，目前主要从事SAR图像杂波建模与变化检测方面的研究工作。广义Gamma模型及自适应KI阈值分割的SAR图像变化检测高丛珊12张红1王超1吴樊11．中国科学院对地观测与数字地球科学中心，北京100862．中国科学院研究生院，北京100049摘要：基于SAR图像的杂波统计特性，利用广义Gamma模型对降噪配准后的SAR图像统计特征进行拟合，获取了辐射值与局部纹理等特征信息；采用信息论中交叉熵的概念，量化不同时相SAR图像统计特征间的差异程度；利用KS与KL检验相结合，自动选取对差异图拟合情况最好的模型，从而实现基于该模型的KI阈值分割。通过对天津市北辰区以南地区的两幅Radarsat图像，以及北京市顺义区的两幅ASAR图像的实验表明，所提出的方法不仅有效地避免了水面波纹变化所产生的大量虚警，并能有效地检测出传统方法所不能识别的，区域内均值不变，仅纹理发生变化的情况。关键词：广义Gamma，K-L散度，KI法则，SAR，变化检测1.引言随分辨率的不断提高，SAR图像越来越多地应用到变化检测工作中。但SAR成像机理的特殊性和数据处理的复杂性，给变化检测带来了诸多困难。尤其是其固有的斑噪、数据获取侧视角度等因素，严重影响了检测结果。因此，该领域自20世纪90年代起，针对SAR图像变化检测主要4四个环节（图1），开展了一系列的研究。其中，以差异图提取与阈值分割的环节最为核心。图像预处理差异图提取阈值分割图像融合图1SAR图像变化检测主要流程差异图提取方面，考虑到乘性斑噪的影响，Rignot（1993）指出比值法较差值法更适用，且对数变换能有效压缩比值图像的变化范围，均值对数比（MRD）因其简单高效的特点，得到了广泛应用。但其对区域内均值不变、纹理改变的情况不够敏感，Carvalho（2001）提出了基于小波变换的差异图提取法；Inglada（2003）引入了信息论中交叉熵（KL）的概念，以计算多时相SAR图像概率密度函数（PDF）间的距离；Mercier（2007）用边沿级数逼近对KL计算进行了优化。目前，仅高斯、皮尔森、瑞利等模型被应用到差异图提取工作中，而实际的地面情况往往十分复杂，简单的模型已不能满足需求，其优化问题成为了该领域的热点。阈值分割方面，传统方法主要分为：最大类间方差法（OSTU）、最大熵法（THC）、最小误差法（KI）等。其中，以OSTU法效率最高，但实验证明该方法在图像目标与背景的面积比很小时失效（Kittler&Illingworth,1986），变化检测方面，变化区域的面积往往有限，该方法不适用；THC法虽可通过圈定阈值选取范围，避免OSTU中目标与背景比例过小的问题，但范围的设定过分依赖先验知识，使结果易随选择的熵准则不同产生很大差异，故普遍性不广；KI法则是Bayes最小错误概率理论的延伸，其代价函数与决策函数随所选拟合模型的不同而变化。其中以Moser（2006）提出的基于广义高斯模型的KI法最为著名。本文将广义Gamma分布模型与信息论中KL散度的概念结合到一起，有效地量化了二时相SAR图像的统计特征差异；并采用了一种KS与KL结合的评价方法，实现对差异图最优拟合函数的自动选取，从而改进了传统KI阈值分割方法。通过对天津市北辰区以南的Radarsat-1和Radarsat-2图像，及北京顺义区的两幅Envisat-ASAR图像进行实验，验证了基于广义Gamma模型的KL差异图提取方法能较准确地量化地表辐射值与局部纹理等特征的变化程度，且通过本文改进的KI阈值分割步骤可获得较之传统变化检测方法，更为理想的效果。2.算法流程提出了一种基于广义Gamma分布模型的KL差异图提取及自适应KI阈值分割的新方法（图2）。2.1广义Gamma杂波统计模型SAR图像固有的斑噪，使传统的差值、比值等方法检测效果不理想。通过求取SAR图像一定区域内的统计分布模型，来削弱其影响；使用两个分布函数之间的差异来表征两幅图像间的差异。SAR图像统计建模技术按其建模过程可粗略分为参量模型和非参量模型（Oliver,2004）。非参量模型建模灵活，适合对未知的复杂概率密度函数进行参数估计，但需大量的样本数据，从而带来了操作复杂、计算耗时等问题。在实际应用中，参量模型更为适用。基于SAR图像乘性斑噪的特点，大多数现有的、广泛应用的统计分布模型，都是由乘积模型发展而来的（图3）。由于SAR图像斑点噪声分量的统计分布是确定的，其统计模型只需考虑表征地物真实RCS起伏的统计分布即可。本文选用的广义Gamma分布（式1）是一个3参数统计分布模型，其中，a是形状参数；b是尺度参数；c是能量参数。形状参数a与能量参数c决定了杂波分布的特性。图像配准&自适应滤波图像配准&自适应滤波原始图像X1(时相T1)原始图像X2(时相T2)利用广义Gamma分布函数拟合原始图像对差异图进行基于选取模型的KI阈值分割将阈值分割结果掩膜与原图进行融合并标注变化信息预处理差异图提取差异图自动阈值分割输出输出变化检测结果图通过计算KL距离生成差异图利用各种经典模型拟合差异图并通过KS&KL检验选取最佳拟合模型利用各种经典模型拟合差异图并通过KS&KL检验选取最佳拟合模型图2算法主要流程采用广义Gamma分布，可以比Rayleigh和Gamma分布更精确地描述拖尾严重的情况，而且可以利用广义Gamma函数随不同的参数估计结果，近似拟合多种不同分布的特性，来对同一幅图像中各子区域具有不同特征的情况，进行拟合。如图4，各种标准的杂波模型，如指数模型（a=1;c=1）、瑞利模型（a=1;c=2）、对数正态模型（a-0,∞）、韦布模型（a=1）、伽玛模型（c=1）等，都是广义Gamma模型的特例。1(;,,)exp()accGpdfcxxfxabcbabb（1）相干斑模型乘积模型斑点噪声模型地物真实RCS分量图像统计模型中心极限定理地物RCS恒定的理想情况推广地物RCS恒定或起伏的情况斑点噪声分量的统计分布RCS的统计分布分解图3由乘积模型发展起来的SAR图像统计分布模型（匡纲要等，2007）图4广义Gamma拟合的经典分布100exp(/)()()(;,,)kGacckkcxxkcadxbabbaxafxabcdxEx（2）最大似然法（ML）是无先验知识情况下的最优估计方法。根据目前的研究现状，无论是利用梯度上升，还是2Pearson距离渐近等算法（Sharifi&Leon-Garcia,1995），其计算量都过于庞大，所以对大尺寸图像进行拟合时，通常采用矩估计的方法快速有效地获取参数估计值。本文采用Gomes（2008）提出的改进的广义Gamma分布参数估计方法，利用k阶矩函数（式2），来获取a、t、c三个参数的估计值。2.2基于K-L距离的差异图提取选用上述广义Gamma杂波模型的概率密度函数对各像元邻域进行拟合，再利用Kullback-Leibler散度，即交叉熵，刻画两个子集之间的相似程度（式3），即通过判断其统计特性的变化，反映两个时相图像区域统计的变化程度（Mercier&Derrode,2004）。()()log()()XXYfxKYXfxdxfx（3）式中()Xfx与()Yfx分别为经矩估计获取参数后，得到的T1与T2两个时相SAR图像中对应像元邻域内的广义Gamma模型。为了解决交叉熵的不对称性，可以采用Johnson-Sinanovic（2001）提出的规则化处理方法（式4）或者Jensen-Shannon（2000）提出的Jeffrey平衡法（式5）解决。(()())/2KLDIVKXYKYX（4）()(()())/2(()())/2QXYfxfxfxKLDIVKYQKXQ（5）2.3自适应KI阈值分割差异图提取后，为在原图准确定位变化区域，需提取变化与未变化类的二值化掩模，经典算法是通过确定变化阈值实现的。本文提出了一种自适应的KI阈值分割方法。KI法则（式7）是从Bayes最小错误概率准则（式6）演化而来。用决策函数()J取代错误概率函数errpz；概率密度直方图()hz则代替先验概率()jP；对数似然概率(,)jpz由代价函数(,)cz表达（式8）。10()|,,j=0,1ZerrjjzpzPpz（6）10()()(,)ZzhzczJ（7）01Z-1,z=1,2,,(,),z=+1,+2,,2ln(),2ln(),zzzppc（8）其中变化类1与未变化类0的后验概率可根据贝叶斯理论，用对应的先验概率和似然概率求取，经简化可得式9。最佳阈值的选取问题即转化为求取决策函数()J最小值的问题，进而获取对应的。0000111110lnln1lnlnZˆˆˆ,;(,)ˆˆˆ,Z-1ˆ()(),(),()()(,)argmin{():=0,1,2,,}iiZzzRzzPpmzPpmPhzzzhzczTcJJ（9）对不同地物，或使用不同拟合模型获取的KL差异图，结果可能近似于不同的分布模型，因此，仅使用经典模型进行拟合，效果不理想。因此，本文引入广义Gauss、G0，及Lognormal和Rayleigh联合分布等模型，同时对差异图进行拟合，并使用KS与KL检验相结合的方法评估拟合优度，自动选取最优拟合函数，并进行基于该模型的KI阈值分割。3.实验与结果3.1实验数据本文使用天津市北辰区以南（图5）、北京市顺义区（图6）两组SAR图像，成像信息如表1。表1实验数据地区天津市北辰区以南天津市北辰区以南北京市顺义区北京市顺义区成像日期20011114200810222004060120050517传感器RADARSAT-1RADARSAT-2ASARASAR极化方式HHHHHHHH分辨率（m）10×810×830×3030×30图像尺寸(pixel)1770×26501770×26501024×10241024×1024（*本文中所有实验选取邻域窗口大小都为7×7）图5天津市北辰区以南(a)2001-11-14(b)2008-10-22图6北京顺义国际机场周边(a)2004-06-01(b)2005-05-173.2差异图提取方法比较均值对数比检测器（mean-ratiodetector,MRD）是目前最常用的多时相SAR图像差异图提取方法。但Ulaby（1986）曾指出，经典的SAR图像纹理分布多数是基于零均值乘性模型的假设。因此，对SAR图像，窗口取得较大时，易发生邻域内均值不变，仅纹理发生变化的情况，MRD方法将无法检测，发生漏检；同样，若窗口过小，对纹理未变化，却因噪声或外部条件变化，使窗口内像元整体变亮或变暗的情况，MRD方法又容易产生虚警。本文利用广义Gamma分布随不同的参数估计结果，可近似拟合多种分布的特点，来描述SAR图像杂波统计特性。采用信息论中用于计算拟合PDF与统计直方图间差异大小的KL散度概念，表征两幅SAR图像的变化程度，从而生成差异图。图7MRD差异图图8G-GammaKL差异图如图7，利用MRD方法对2001年与2008年的两幅天津市的Radarsat图像进行检测，所获取的差异图水域部分产生了大量虚警；利用基于广义Gamma分布的KL距离获取差异图，如图8，可有效地避免这一问题。为了进一步验证基于广义Gamma分布的KL差异图提取方法的优越性，本文从图9中选取了对应的三块样区，A为西青区鸭淀水库、B为津南水库、C为北大港湿地。已知样区大小为1186048个像元，其中实际变化像元有174608个，未变化像元有10