基于单边拟合策略的合成孔径雷达图像变化检测毕业设计.doc

致 谢35 第一章 绪论1.1 SAR成像基本原理及图像特点雷达主要是利用电磁波的散射进行目标探测和测距，具有全天候和强透射性等优点，使其在各个领域都发挥着重要作用。传统雷达主要解决目标有无的问题，获得回波脉冲及强弱，并不产生图像。后来人们利用多普勒锐化技术，实现高分辨率的雷达成像。SAR是在雷达基础上发展起来的一种探测技术，它是一种利用雷达与目标的相对运动把尺寸较小的真实天线孔径用数据处理的方法合成较大的等效天线孔径的雷达。为了获得二维图像，在距离上采用线性调频信号技术，可以从其回波中将较近的目标分离出来，从而达到高分辨率;在方向上，通过平台的运动构建一个虚拟的大孔径天线，可以获得对另一维的高分辨率1。合成孔径技术的基本思想是对在平台前进方向的不同位置上所接收的包含相位信息进行记录和处理，从而得到比实际天线更长的假设天线进行观测的结果2。SAR获取图像的几何原理如图1.1所示。图1.1 SAR成像原理示意图成像处理一般包括三个步骤：1、距离向脉冲压缩，对一个合成孔径内的每个回波脉冲均进行压缩；2、距离徙动矫正，把属于同一点目标的脉冲校正到等距离线上、3、方位向脉冲压缩。难点在于方位向，需要利用回波信号的相位提取多普勒频率。平台轨道和姿态的不稳定及测量误差、雷达系统的不稳定等，均会引起信号的幅度和相位误差，从而影响匹配滤波的处理精度，很难实现理论的分辨率（D/2）。方位压缩的滤波匹配需要估计出两个参数：多普勒质心和多普勒频率。常用的算法有：距离-多普勒（Range-Doppler,R-D）算法。 Chirp-Scaling(C-S)算法SAR主要采用侧视雷达，以有源主动方式工作，其特殊的成像方式，形成了图像的独有特征，包括：辐射特征、噪声特征、目标集合特征、统计特征、纹理特征等。雷达图像的分辨率有距离向和方位向两种。雷达图像的分辨率和像元大小是两个概念，雷达图像中每个像元大小通常低于分辨率。SAR图像的灰度统计特征：雷达图像上色调的变化，主要取决于目标物的后向散射截面。每一个接受到的回波被转换成电信号，并以某一特定的灰度色调记录在光学胶片上或转换成一个具有特定值、由于表示亮度的数字化像元。坡度的变化，复介电常数（含水量）、表面的粗糙度 是影响雷达图像色调的三个主要因素。其中表面的粗糙度在决定雷达图像的灰度（也就是回波强度）上起着决定性的作用。SAR图像与传统的光学图像不同，它是通过天线发射有规律的无线电信号，并接收从地面返回的回波信号然后采取用成像处理的方法将回波信号转化为实际图片。雷达成像是依靠对回波信号提取目标的后向散射系数，反映所探测区域的微波特性，光学成像依靠的是普通的光学反射1。相对于与光学图像，SAR图像具有许多重要特性：1 全天候、全天时，受风雨天气影响非常小，不受太阳照射影响，白天黑夜均可工作2 对电介质（水分量、生物量、冰雪等）、地表粗糙度（海洋风速等敏感）3 精确测量距离4 对人工目标以及目标结构较敏感5 强透射性，能穿透遮盖物，比如植被、树叶6 分辨率不受雷达波长、载体飞行高度、雷达作用距离无关1.2 SAR技术发展历程雷达卫星是载有合成孔径雷达（SAR）的对地观测遥感卫星的统称。尽管迄今为止，已在一些发射的卫星上携有SAR，如Seasat SAR, Almaz SAR, JERS-1 SAR, ERS-1/2 SAR, 与它们搭载在同一遥感平台上还装载着其他传感器。而 1995年11月发射的加拿大雷达卫星（Radarsat）则是一个兼顾商用及科学试验用途的雷达系统，其主要探测目标为海冰， 同时还考虑到陆地成像，以便应用于农业、地质等领域。该系统有5种波束工作模式，即：（1）标准波束模式，入射角2049 ，成像宽度100公里，距离及方位分辨率为25米x28米；（2）宽辐射波束，入射角2040 ，成像宽度及空间分辨率分别为150公里和28米x35米；（3）高分辨率波束，三种参数依此为3748，45公里及10米x10米；（4）扫描雷达波束，该模式具有对全球快速成像能力，成像宽度大（300公里或500公里），分辨率较低（50米x50米或100米x100米），入射角为2049；（5）试验波束，该模式最大特点为入射角大，且变化幅度小4959，成像宽度及分辨率分别为75公里及28米x 30米。合成孔径的概念产生于20世纪50年代初期。1951年6月，美国Goodyear航空公司的Carl Wiley首先提出利用频率分析的方法来改善雷达角分辨率，并称其为多普勒波束锐化。1952年C.W.Shervin提出了全聚焦陈列概念以及运动补偿概念。这些思想促进了X波段相干雷达的研制。1957年密歇根大学雷达与光学实验室研制的SAR系统获取了第一张全聚焦的SAR图像。50年代末，美国成功研制第一批可供军事侦察用的机载高分辨力合成孔径雷达。1967年Greeberg提出在卫星上安装SAR的设想。从此，人们开始了一航天飞机、卫星等为载体的空载SAR研究。60年代中期，随着遥感技术的发展，合成孔径雷达技术从军用推广到民用方面，并发挥了重要作用。1978年，NASA发射了全球第一颗装载了合成孔径雷达的人造卫星SEASAT,此卫星标志着SAR进入太空安装SAR对地观测的新时代，星载SAR由实验阶段转向实际应用的阶段。1981年，美国“哥伦比亚”号航天飞机搭载SAR-A 顺利升空以及1984年进行了“挑战者”航天飞机搭载SAR-B的实验。1987年前苏联发射了第一部长期运行的空间合成孔径雷达，名为ALMAZ-I。1988 年12 月美国用“阿特兰蒂斯”号航天飞机将世界上第一颗高分辨率雷达成像卫星“长曲棍球”(LACROSSE)SAR卫星送入太空,其分辨率已达1m,在海湾战争中发挥了重要作用。欧空局于1991 年7 月和1995 年4 月分别发射成功ERS-1和ERS-2卫星。日本于1992年发射了JERS-1空间合成空间雷达，工作于L波段。1994年9月，美国航天局、德国空间局和意大利空间局共同进行了航天飞机成像雷达飞行任务。1995年，加拿大发射了第一颗资源勘测卫星Padarsat-1，该系统为商业和科学研究提供了全球的冰情、海洋和地球资源等数据。SAR成像技术随着时代不断前进发展，且功能实现更加完善，包括频率方面由单频段到多频段、成像空间由单极化到多极化，视角角度方面升级到多视角，性能的大大升级扩大了对全地球的覆盖率，提供的全日时的合成孔径雷达探测数据信息准确度更加准确，SAR逐渐向小型化、SAR星座组网及干涉测量等方向发展。目前合成空间雷达分辨率已达到0.1m数量级，纵观国外SAR的发展历程，可以看出随着科学技术的不断进步，SAR的水平和功能也在不断提高。可以相信，科学家们将不断地挖掘SAR的技术潜力，为人类的需求服务。我国SAR研究开始于20世纪70年代后期，起步较晚。经过几代人的不懈努力，许多技术已经走在世界前列23。1.3 SAR图像变化检测介绍SAR图像变化检测是指对不同时段同一地区的遥感图像，经过一定的分析和运算得出地面的变化区域、特征以及过程，涉及变化的类型、分布状况与变化量。SAR图像变化检测的主要任务是检测局部文理和辐射值的变化。目前SAR图像变化检测思想主要分为两大类，一种是非监督的变化检测思路。常用非监督的变化检测方法常用的例如图像比值法4、图像差值法。此类方法优点在于简单、直接。如图1.2所示差，异图分类法的关键在于差异图的构造和阈值的选取。 时相1图像时相2的图像图 像预处理图 像预处理比 较处 理变化检测结果图差 异图 像差异图像分析 图1.2 非监督的变化检测方法另一种主流思想是基于分类的变化检测方法思路。如图1.3所示首先对不同时段的图像分别进行分类处理，然后根据分类结果比较，得出类别的不同并获得变化的区域。后分类比较法的关键在于多时相图像的准确分割和分割结果与差异图变化信息的提取与反馈。时相1的图像时相2的图像图像预处理图像预处理自动分类处理变化检测结果图分类图分类图比 较处 理自动分类处理图1.3 基于分类的变化检测方法时相1的图像时相2的图像图像预处理图像预处理自动分类处理变化检测结果图分类图分类图比较处理自动分类处理由于星载SAR遥感术后局具有全天候、全天时的独特优势，同事具有固定的重访周期，所以它比光学遥感更适合于变化检测。针对SAR图像的特点，目前比较通用的方法主要有差值法、相关系数法、直方图比较法、图像熵方法来对SAR图像来进行变化检测2。归结起来，近些年逐步形成的SAR图像变化检测方法可以分为：(1)基于简单代数运算的变化检测，经典的变化检测方法包括图像差值法、图像比值法、对数比值法；(2)基于图像变换的变换检测，经典的变化检测方法包括主成分分析、变化向量分析法、相关分析法图像变换法；(3)基于图像分类的变化检测方法。1.4 SAR图像变化检测发展现状及应用国外在应用SAR遥感进行变化检测方面起步较早，进行了一系列研究。Cyril采用模糊隐马尔可夫链方法进行SAR变化检测5，通过计算两幅配准数据的对数比值，并引入模糊隐马尔可夫链，最大限度地提取相关性进行变化检测。Ramin6提出来通过模拟部分极化电磁场引入向量和观察量而建立变化空间，从而对相关或变化特征分解，更深入理解图像变化的本质。Andreas7在小波变换和位置上衣脊状特征重建原始图像，在曲波空间进行变化检测。Hachicha和Chaabane8应用Dezert-Smarandache理论进行了变化检测的探索。Eric采用ERS-1数据对SAR变化检测进行多种方法的初步研究和对比9。Onana针对ERS-SAR多时相图像的特点，采用了数据融合的方法解决城镇地区的变化检测10。Patricia和Kenneth分别对SAR相关系数的变化检测方法进行了深入研究11,12。Francesca将小波变换的多尺度方法引入到SAR图像变化检测中13。Jordi将聚合交叉熵计算方法引入到SAR图像变化检测中，取得了初步的结果13,14。在国内，由于SAR数据源不足的问题，开展SAR变化检测的研究相对数少，近几年来的主要成果有：基于二阶灰度统计特性的SAR图像变化检测15，该方法着重研究了统计特征中的二阶灰度特征如二维直方图、自相关、能量、协方差、惯性矩、熵、绝对值等，先对多幅SAR图像分别计算灰度特征，然后用差值法对不同时相图像内相应的特征进行变化检测，最后用数学期望最大算法获得阈值。基于分形维数的SAR图像变化检测16，引入盒子维计算法，该方法最大的优点就是对SAR图像固有的斑点早说不敏感。也有学者引入马尔科夫随机场理论，研究了一种顾及空间领域关系的多时相SAR图像变化检测方法，该方法能有效提高漏检率，实现变化特征不明显区域的有效提取17。对于现在SAR图像变化检测方法大多根据SAR幅度的不足，有学者针对SAR图像的极化特征及其分布模型构建极化似然比检测模型18，实验证明该方法能有效地区分地物变化情况。针对传统分类比较法存在的分类误差累积问题，根据不同时相图像的不变信息的相关性，有学者提出了一种基于两时相图像联合分类的SAR图像变化检测方法19，该方法以灰度值作为输入信息，通过相似度计算得到两时相图像对应位置像素的灰度相似度，然后求解全局相似度阈值，并用K均值对两幅时相图像进行联合分类，最后通过类别比较获得变化检测结果。由于SAR图像变化检测技术在国内起步较晚，加之SAR数据源不足，因此在变化检测的方法和技术上没有根本的突破，误捡和漏检现象比较严重，大大影响了检测结果和精度。如何使变化检测技术更具稳健性和通用性，如何更好地利用图像中感兴趣目标的结构信息，降低图像中噪声和匹配误差带来的影响是未来工作的重点。变化检测具有广泛的应用。例如，监测森林和植被变化；监测农作物的生长变化、土地覆盖利用变化；监测和评估各种灾害发生前后变化；军事目标区域的动态监视、打击效果评估等。1.5 本论文的主要工作本论文主要在前人的研究的SAR图像变化检测方法基础上，设计一种基于单边拟合思想，单边拟合差异直方图，自动确定阈值的SAR图像变化检测方法。本论文的工作分为四章，各章内容如下：第一章介绍了合成孔径雷达成像的基本原理和SAR图像的特点，简述了SAR技术发展历程，包括国内外发展情况及各国成就，简单概述了SAR图像变化检测，包括变化检测方法、过程和在各个方面的应用。在此基础上，了解了研究目标所需要的各个细节及知识，并为论文工作打下良好基石。第二章主要介绍合成孔径雷达图像差异图的构造。本论文研究方法是一种无监督阈值方法，此方法首先要对两时相遥感图像进行比较，获取差异图，在此基础上，设计变化检测方法。差异图的构造方法有多种，本章主要介绍各种差异图构造方法。第三章主要就该SAR图像变化检测进行详细描述。在获得差异直方图后，设计相关的模型，结合相关参数估计策略求得变化类和非变化类的概率统计分布，利用贝叶斯决策理论自动地确定阈值。本方法结合了EM算法等相关算法，并作详细描述。第四章总结本论文的工作，并指出不足之处与优点，以及以后的研究方向。第二章 基于差异图的构造SAR图像变化检测方法SAR图像变化检测的过程可以简单概述为3步：1.构造差异图，2.对差异图进行分析与处理，3.生成变化检测结果图并分析。构造差异图作为第一步，是开端，起着至关重要的作用。对于差异图的构造方法，众多学者纷纷研究，并提出许多算法。目前广泛运用的算法主要为差值法和比值法。2.1 变化检测技术存在的问题变化检测技术在图像处理领域中是一个活跃的研究热点，影响其检测效果的因素主要有以下几点：u 图像配准许多遥感图像的应用都需要求对图像进行配准，列如变化检测、图像融合、图像镶嵌、基于两幅图像提取DEM等。配准是将两幅遥感图像上的同名点准确配准的过程，优势图像之间是相对配准的，此时像素地标可能是未知的。有时配准过程中将会出现大量虚警现象，这造成的原因可能是由于没有经良好的精准配准处理而导致的图像错位造成的。目前，在导致图像变化检测的精度出现错误率的原因中，对于图像配准的精准度方面的研究成果与分析，还是比较少的。Dai 和 Khorran21关于此方面做了一个相关研究：采用一个多频谱图像，对其连续进行自我变化处理，然后采用基于差值的图像变化检测处理，并同时分析其与图像配准的精确度的相关性。他们得出的结论是：好的变化检测结果需要高精度的配准算法，如1/5个像素的配准精度获得的变化检测的错误率小于10%，也就是说配准精度越高得到的变化检测效果越好。u 阈值的选择为判断区域是否发生了变化现象，许多变化检测方法会通过构造阈值的方法，来区分变化区域和非变化区域以得出。目前构造阈值的主要方法有两种：第一种为交互式方法，第二种是通过采用一些统计检测量来自动地获取阈值的方法。前者是通过人为自己调节阈值的选取直到取值得到满意或者进行检测时采取经验值处理；后者为非监督变化检测阈值获取方法，这是变化检测技术未来研究发展的一个重要方向。虽然目前阈值的获取方法有多种，但是最优获取阈值依然是目前在变化检测方法研究领域中一个重点，也是需要攻克的难点。u 辐射校正辐射值的改变可能是由于感兴趣的变化导致的，也有可能是其他一些随机原因造成的。后者包括的因素主要有大气状况、辐射的角度大小以及地区土壤的湿度等因素。在进行变化检测时，要确定一个条件，那就是前者的造成的改变要比后者造成变化要大。如果后者造成的变化较大，走就必须要求进行辐射的校正。在进行辐射校正时，一般会对多时相遥感图像的非变化区域进行辐射校正，但在目前实际的应用中，大家会利用一幅遥感图像作为参考图来实现辐射校正。在基于合成孔径雷达图像的图像变化检测中，我们可以将两幅图的灰度值进行匹配，来缩小图像之间的灰度差值。u 性能评估目前对变化检测算法的性能评价并没有一个统一的衡量标准，根据不同的应用需求大体分为定性评价和定量评价两种。定性评价也称为视觉评价，大家通过计算机卷帘或闪烁显示技术实现，适合于人眼的直接观察。定量评价在实际应用中具有很大难度，究其原因是较难获取真实的地物变化信息。其中多余变化检测进行效果评价时 对地表数据获取真实情况是其中一个难度。在不同时间或不同情况下，对于一幅图像的评价，一个人的判断可能获得不同的结果。u 感兴趣变化在具体应用中的定义对于不同的情况，或不同的实际应用，每个人都有可能得到对变化区域的不同定义，变化的区域可能是由非变化的区域或是变化的区域定义得来。定义变化的区域不同，对不同的SAR图像变化监测可能构成不同的影响。2.2 SAR图像的预处理在进行对SAR 图像的一些操作时，如果直接用原图操作，可能对后续步骤造成不好的影响，所以在进行变化检测前，需要对图像进行相应的预处理措施，主要是因为合成孔径雷达进行多时相成像时可能是在不同条件下得到，而且影响合成孔径雷达成像效果的因素包括受一定的斑点、噪声等影响，这些不良因素对于进行下一步的SAR图像变化检测算法构成大量的困难以及会降低图像的效果质量。造成SAR 成像出现的斑点噪声等不良效果是由于SAR使用采用一些相干信号而导致的，这是成像系统自身所自带的。SAR 图像的斑点噪声不良现象表现为图像上的亮度起伏不定，这会让人感觉图像很乱，而且探测目标的边缘表现不够清晰。因此，降低图像固有斑点和噪声的措施在图像变化检测中是非常重要的。现在的相干斑的抑制技术可以大致分为两种，一种为多视技术，另一种为对 SAR 图像空间的滤波技术。前者技术在以前的 SAR 处理采用率较高，利用空间的分辨率来抑制相干斑是此方法的重要特点，此方法需要将空间的分辨率减小到之前的1/ N，且相干斑抑制变为之前的。此方法主要基于所有的像素点的像素值完全相等的思想，缺点在于忽略SAR图像边缘结构的信息，而后者技术在SAR图像处理中则是基于在尽量保证以不减小空间分辨率作为前提下，加强抑制相干斑。Lee 滤波、Frost 滤波、Kuan 滤波是比较经典的几种SAR图像滤波方法。除滤波外，对图像进行高精度配准，在SAR图像变化检测中也是非常重的。在实际应用中，许多变化检测方法都要求高精度的图像配准，否则较低的配准精度对导致出现大量虚警现象。图像配准是图像变化检测中一步，但是不可忽略，应该说这也是人们需要研究的一个重要领域，因为配准的精度高低会影响SAR图像变化检测的效果。在图像配准方法中，基于图像特征点来进行图像配准是最为广泛的方法。但是这类方法的缺点往往无法很好地确定特征点的位置，主要由SAR 图像中斑点噪声影响严重造成。因此在 SAR 图像的配准过程中，我们往往会通过人为手动选取特征点的方法来对图像进行精度配准。预处理有两种形式，有的预处理方法需针对两幅图像进行，而某些方法以其中一幅图像为基准，只对另外一幅图像进行处理。SAR 图像的灰度匹配主要是指通过对SAR原始图像的灰度进行变换，消除由传感器扰动，成像距离差异等随机原由所带来的影响。在光学图像中一般采取的方法主要有强度归一化、同态滤波和图像回归法。而在SAR图像变化检测中，主要适用的方法为第一种和第三种。u 强度归一化在 1973 年，Ulstad22提出这种方法用来消除数字图像中差异现象。对成像距离进行等变化操作会对灰度直方图进行平移、伸缩，结果会使多时相两幅图像的均值和方差变得不相同，而通过此强度归一化方法可使得两幅图像会具有相同的均值以及相同的方差。设图像对为和,和分别表示两幅图像中对应同一点处的幅度值，利用和进行强度归一化得： (2-1)其中,分别为和的方差和均值。u 图像回归法此方法假设将检测图像的像素值映射为参考图像与像素值相对应的一个线性函数，之后可以通过最小二乘方法来估计此线性函数。这种方法主要对应于不同时间段内图像的像素间的方差和均值存在差别的时候进行处理。具体公式如下： (2-2) (2-3)式中，是线性函数中参数且待测定。一般的处理方法是选择一些已经确定像素值没有变化的点，再利用这些点来估计参数值，。2.3 基于比值法的SAR图像变化检测算法2.3.1 算法介绍首先对两幅变化前后的SAR 图像分别进行滤波处理，消除噪声，尽可能的减少噪点对变化检测结果产生影响。然后利用图像比值法生成初步的新像素矩阵。假设是第一幅变化前SAR 图像对应像素点的灰度，是第二幅变化后SAR图像对应像素点的灰度。对应两幅图像，变化的部分像素差异波动较大，所以比值/与1距离就比较远；而对应不变的点，像素虽然不至于完全一致，但是比值应该接近1。这个特性就为此方法的实现奠定了良好的基础。具体公式为 （2-4）为了克服噪声，也可以采用一定大小的窗口进行平均后再进行比值运算，即 （2-5）比值操作后所得到差异图的统计分布如下面公式所示： （2-6)其中,表示在两个时刻同区域的均值图像强度，N表示SAR图像数据的等效数据。然后对上一步得到的像素矩阵进行处理，对比值如进行取反操作，再乘以某一恰当的系数将像素值的差距拉大，便可得到一幅SAR图像变化图像。最后将实验所得的结果图跟已知参考SAR图像进行比对，获取变化区域及变化特征，并统计漏检和误检的像素点数量。2.3.2 算法步骤步骤1、首先对两时相同区域SAR图像进行几何校正和辐射校正步骤2、通过基于比值算子产生差异图步骤3、产生变化结果图并分析2.4 基于差值法的SAR图像变化检测算法2.4.1 算法介绍图像差值法的基本原理是首先计算不同时相图像对应像素灰度值或纹理特征值的差值，生成差值图像，接下来对差值图像选择合适的阈值找出差异较大的部分，以表明此期间地区陆地表面的变化。差值法是简单而直接的变化检测方法，首先对两幅变化前后的SAR 图像分别进行滤波处理，消除噪声，尽可能的减少噪点对变化检测结果产生影响。将两幅图像的对应像元相减，就可以得到一幅差分SAR图像，即 （2-7)为了克服噪声干扰以及提高算法有效特性，常常采取一定大小的窗口进行平均，应用均值进行差分处理，即 （2-8)式中：I1,I2 -时相1,2的SAR图像上的像素；n-窗口的像元素。阈值对于图像差值法非常关键，需要凭借经验的累积和多次试验从而选取合适的阈值。图像差值法在理论上的最直观并且非常实用，是目前应用最广泛的变化检测方法。2.4.2 算法步骤步骤1、首先对两时相同区域SAR图像进行几何校正和辐射校正步骤2、通过基于差值算子产生差异图步骤3、产生变化结果图并分析2.5 实验结果差异图的构造实验，本论文取加拿大Ottawa地区的1997年5月和1997年8月Radar-sat影像，为两时相的SAR图像。首先对两幅图像预处理，再采用3种差异图算法相对应地构造差异图。为了比较各差异图的优越，通过构造变化检测结果图，并进行对比，比较各算法的差异。 实验原图 (a) (b) (c) 图2.1 Ottawa SAR图像前后变化图像及变化检测参考图图2.1中（a）、（b）为SAR图像原始图像，分别为加拿大Ottawa地区的1997年5月和1997年8月Radar-sat影像。两幅图像的大小均为290350，灰度级为256。其中，前后时段图像的变化信息主要是由于夏季雨季来临，洪水淹没部分陆地区域所致。图像（c）为变化检测参考图，白色区域表示变化区域。 差异图 (a) (b) (c) 图2.2 Ottawa经3种算子构造的差异图图2.2表示Ottawa两时相SAR图像经各种算子构造的差异图。图（a）差异图由差值法构造，图（b）由比值法构造，图（c）由对数比值法构造。 检测结果图 (a) （b） (c) 图2.3 Ottawa模糊C均值聚类法后的变化检测结果图构造差异图后，利用模糊C均值聚类法23，将差异图像素点分成两类，构造变化检测结果图。图2.3中3幅图分别为差值法、比值法和对数比值法的变化检测结果图。结果分析：通过比较3种方法构造的变化检测结果图，可以看出，图像差值法构造的检测结果图漏检数较多，而图像比值法构造的检测结果图虚检数较多，相对比对数比值法构造的变化检测结果图与参考结果图更接近。图像差值法算法较为简单，易于理解，且最直观，但图像变化区域漏检较多。对数比值法相对比值法，能把图像固有的乘性相干斑噪声模型转化为加性噪声模型，大大提高了比值法的精确度，对于实际应用有很大意义。第三章 基于单边拟合策略的SAR图像变化检测本论文设计的SAR图像变化检测方法，与基于差值法和基于比值法的SAR图像变化检测不同，它是在差异直方图的基础上，构建变化类和非变化类的概率分布图，并拟合直方图曲线，且只限于单边拟合，使得当拟合曲线与直方图曲线在单边区域达到平行。所以要结合相关的模型参数估计策略来获得变化类和未变化类的概率统计分布。得到概率分布函数后，需使用算法来自动确定最终的阈值，此时需要应用贝叶斯决策理论，并利用此阈值来构造变化检测结果图。本章的工作需要掌握相关的概率模型知识，以及如何运用贝叶斯决策理论。3.1单边拟合原理假设（=0，l，.-1)表示为经处理之后的差异直方图。我们可以将h(x)看作这个差异图的概率统计分布图。确定差异图的阈值实际上就是区分差异直方图的未变化类灰度值与变化类灰度值。差异直方图形式可以表现为三种，分别为理想的差异直方图、无混叠的非理想差异直方图以及混叠的非理想差异直方图。图3.1所示理想差异图直方图曲线，这种情况主要基于没有任何的几何与辐射校正误差、斑点噪声的影响。图3.2与图3.3表示为非理想差异直方图，可以看出这两种差异图相对于理想差异图，表现为差异直方图曲线分散和平移，导致此现象为几何与辐射校正误差、斑点噪声等影响。图3.2表示为，在变化类与非变化类之间有很少的混叠，而图3.3表示为，变化类与非变化类之间有着较为严重的混叠，此情况较为常见，但是此情况的阈值相对于前者较难以确定。较为常见的处理方法都是先构造直方图中变化类与非变化类别的概率统计模型。一般的概率模型很难与差异直方图进行拟合，主要是由于两方面：一是差异图直方图情况复杂多样性，二是现在所构造的概率统计模型具有局限性。根据此情况，不考虑完整拟合，而只考虑直方图中变化曲线和非变化曲线进行拟合。从下面3幅图像中可以看出，非变化类以Mu点分散，而变化类以Mc点分散，所以阈值的产生一定在Mu值和Mc值之间。在本论文设计的方法中，曲线拟合值针对直方图曲线的一边，而不同是拟合两边，这样可以大大提高直方图曲线拟合的精度。 图3.1 理想变化检测直方图图3.2 非理想变化检测直方图（无混叠）图3.3 非理想变化检测直方图(混叠）本论文设计的单边拟合策略，只将Mu与Mc之间直方图曲线与拟合曲线进行拟合。在进行拟合的时候，只进行对Mu右侧单边非变化曲线的拟合和对Mc左侧单边变化曲线的拟合，其它区域的曲线不进行考虑。通过直方图的信息可以看出，变化区域像素相对于非变化区域的像素较少，这与实SAR图像相符合。由于此情况，使得变化区域的曲线坡度较低，有时很容易被混叠区域覆盖。为了简化算法，本方法不考虑变化区域的单边拟合，而直接求出变化区域的条件概率分布。3.2 基于广义高斯模型的SAR图像变化检测模型3.2.1 广义高斯模型的建立在SAR图像变化检测中，广义高斯模型往往被应用于构造条件概率分布。为了满足单边拟合的条件，则需要在广义高斯模型加入一个单边平移参数，构造出相对应单边拟合策略思想的广义高斯概率模型。假设差异影像(0)上，表示最大灰度值)变化类为，非变化类为。所以单边曲线的条件概率分布函数可描述如下： (3-1) (3-2) (3-3) (3-4)参数,分别是单边广义高斯分布的方差、形状控制、均值和单边平移参数。(3-3)式中为Gamma函数。通过改变添加在广义高斯分布函数的单边平移参数大小，移动概率分布曲线，来实现非变化单边曲线的曲线拟合。从上面式中可以看出，在单边广义高斯模型中还需要进行估计的参数主要包括未知的形状参数与单边平移参数。3.2.2 广义高斯模型的参数估计 参数估计设差异图像的灰度动态变化范围为0, 1，(=0,1,.,-1）是差异图像上任意像元个数占像元总数的百分比，表示灰度值的概率，再设选择的初始化分割阈值为。 广义高斯模型中的形状参数的估计方法很多，考虑到估计过程的简单性与参数估计的精确性，本文使用文献24提出的估计方法估计参数，该方法通过样本的方差和期望就能估计出参数，而且在文献24的实验中也体现了该方法比其他方法更为精确。 平移参数c的估计本论文设计的单边广义高斯模型的思想是先构造非变化区域的广义高斯单边曲线与非变化单边曲线进行曲率拟合，最后通过比较，选择最优的非变化曲线上，并得到此种情况下的概率统计分布。在构造广义高斯单边曲线时，会将广义高斯单边曲线和非变化单边曲线进行比较，当两个曲线似平行时，可以认为此时的曲率拟合是非常理想的，同时我们发现两者曲线之间会存在一定的间隔，结合本论文单边拟合思想，为了得到最佳拟合，就需要移动广义高斯单边曲线使得两者曲线一部分相重合。估计广义高斯模型中单边平移参数可以具体概括为下面两个步骤：1首先确定直方图中非变化单边曲线的灰度值范围。确定区域间步骤如下： (1)首先确定非变化单边曲线的区域初始值，及下限值为等于直方图上的Mn值。 (2)求解非变化单边曲线的区域的上限值。由于直方图上变化曲线与非变化曲线可能产生混叠，从而掩盖重合区域，因此并不容易得出上限。从图3.2和图3.3，我们可以看出，当混叠区域并不明显时，变化类和非变化类曲线呈现双峰形状，而混叠区域非常明显时，非变化曲线相当于变化类曲线的一部分。由于单边曲线的复杂性，我们要根据不同的情况，来确定相对应的非变化曲线的上限值。计算非变化单边曲线上限公式如下： （3-5)其中，e=10-6，T为直方图平滑因子，本章中取T=10。如果满足公式 （3-6)即曲线呈现双峰状则单边拟合区域没有产生较大混叠，否则产生较大混叠。(3)判定非变化区域的范围。如果产生了混叠现象，则非变化单边曲线区间范围为，如果混叠现象没有产生，则非变化单边曲线区间的上限不改变，下限变成，这样增加下限的值，也导致单边曲线拟合区间范围减小，这样有利于提高单边拟合的精度。2计算平移参数。主要由下面步骤完成：(1)初始化=0，将估计的参数的估计值代入式（3-2）、式（3-3）和式（3-4）中，分别得到参数和两值，让后将4个参数、代入式中（3-1）中(2)利用下面公式求解单边平移参数： （3-7)其中，代表向量的均值，参数表示方差，参数表示均值。3.3 基于最小错误率的贝叶斯决策的自动阈值选取法人们希望模式分类问题中尽量减少由于分类出现的错误，为了能达到这样的目的，可以通过利用贝叶斯公式在概率论中应用，然后就可以使得分类错误率变得最小，这种分类方法策略称为基于最小错误率的贝叶斯决策理论14。在SAR图像的变化检测过程中，我们将差异图像上的像素点分成变化类和未变化类两种，此方法的阈值选取就是要找出基于采用错误率最小的分类方法来将变化的像素点和未变化的像素点区分开来。差异图像上的像素点分为变化类和未变化类，分别用和表示，和分别表示为各自分布的先验概率，随机灰度值在和中出现的概率，称为似然概率密度函数。差异图像的直方图可以被看作是概率密度函数的估计。差异图像总概率密度函数是由变化和未变化两类概率密度函数混合组成，其表达式为 （3-8)利用贝叶斯公式有: ，其中 （3-9)条件概率称为后验概率。结合贝叶斯公式，观察灰度值，利用某类先验概率通过公式求出该类后验概率。这样，我们总结最小错误率贝叶斯决策的规则，可以概括如下：如果 ，则将归类于未变化类 ，反之 则将归类于变化类。利用贝叶斯公式还可以得到几种最小错误率贝叶斯决策规则的等价形式：1.如果，则。2.若，则，否则那么最佳阈值值，应当满足下式: （3-10)先求得似然概率密度函数，如果两个表达式己知或者可以通过假设来获得，则它能够确定一个具有最小的错误率的最佳门限，将图像中变化和未变化的像素区分开来。但在实际中，似然概率密度函数是未知的。因此，我们一般的做法是结合相关图像分布模型，也可以利用样本参数的方法来估计概率分布，接着模拟出似然概率密度函数，然后估计模型参数，利用这些参数得到似然概率密度函数，最后利用这些函数求解要得到的阈值。3.3.1 KI 最佳阈值选择算法在1986年， Illingworth和Kilter在文献25中提出的一种阈值选择方法，且基于最小错误率贝叶斯理论，称为KI阈值选择方法。它的主要思想在于为求得图像分割的最佳阈值，当针对服从一定分布模型的图像时，采用一种简单的参数估计方法获得。差异图像中，如何选择阈值，KI算法的运用具体介绍如下。 设通过算法选择的阈值为,则差异图像中大于的部分认为是变化部分，而其余部分为未变化部分。我们定义、用来描述未变化类和变化类的先验概率、均值、方差，而相应的公式如下，其中为灰度值的像素出现的概率： （3-11) （3-12） （3-13） （3-14） （3-15） （3-16)定义一个函数为准则函数，用来描述整个图像上平均的正确分类性能: （3-17)其中为代价函数，它是一种度量，反映分类性能，其表达式为: （3-18)表示为在给定阈值和灰度值时，未变化和变化类的后验概率密度。选择最佳阈值的目的在于要使图像的平均分类性能达到最小值，即: （3-19)根据贝叶斯理论，可以用先验概率和似然概率来表示后验概率密度函数，即 （3-20）准则函数可变为： （3-21）其中表示类别的熵, 。具体表达式为： （3-22）使得达到极小值的值即为最佳阈值。3.3.2 EM 迭代阈值选择算法KI算法的思想是通过使反映整幅图像正确分类的准则函数达到最小，来找到最佳阈值以区分变化和未变化类像素点。相对比，当变化类和未变化类的似然概率密度函数的数学期望达到最大值时，是EM算法选择最佳阈值的判断条件。最大数学期望算法(EM: Expectation Maximization Algorithm)不需要任何先验知识和外来的数据，是对不完整数据问题进行最大似然估计的一种常用算法，只需要从观测数据自身来获得参数的估计值。标准的EM算法可以看做一个迭代的过程，每次迭代均需要两个步骤组成，包括先求期望值和之后期望值最大化，前面一步通过待估计参数的当前值，再结合观测数据中来直接估计概率密度的期望值大小，后面一步通过最大化这一期望值来重新估计参数的量，在迭代过程，两步步骤是依次交替进行的，直至收敛结束迭代。基于高斯模型的差异图像上未变化的像素点和变化类的像素点分别为和，同时相对应的似然概率密度分别为，其表达式表示如下: （3-23），分别表示为未变化的像元类和变化的像元类的均值以及方差。在上述假设条件下，在上述算法的公式中，将有4个参数未知，这可以利用差值图像来估计得到。通过EM算法，计未变化像元类统计分布参数。差值图像的行数和列数表示为和，差值图像中第行第列的像元灰度值表示为，类条件概率的计算公式为(3-22)，全概率的值通过式（3-5）给出。而且变化类的相关参数先验概率,均值和标准差 的计算方法同上。当阈值的计算终止时，其条件为当相邻两次迭代计算的，和的值小于给定的阈值。3.4 算法步骤步骤1：取得两幅不同时间同地域的SAR图像，预处理，采用Lee滤波，差值法构差异图像步骤2：求出差异图像的灰度直方图，确定单边拟合区域，区域的范围步骤3：令初始阈值=Mu,构造单边广义高斯模型，求出阈值Ti下的无变化区域直方图曲线的概率分布函数步骤4:利用准则函数计算阈值下的单边拟合误差,令阈值=+1，重复计算阈值步骤5:令阈值产生区域中拟合误差最小的阈值为，确定无变化区域的直方图概率分布步骤6：利用高斯模型与期望最大化EM参数估计方法，求出变化区域直方图曲线的概率分布函数步骤7：根据最大后验概率方法确定最终阈值，并通过该阈值构造变化检测结果图3.5实验结果l 实验原图 （a） (b) (c)图3.4 Bern SAR图像数据集原始图像及变化检测参考图图3.4中前两幅SAR图像分别为1999年4月和1999年5月通过ERS.2获得的瑞士Bern地区的前后两时相SAR图像组成，时相SAR图像（a）是在拍摄于洪水灾害刚刚发生后的情形，图像中像素有高低，其中较暗的部分为受洪水影响淹没的地方，时相SAR图像(b)是在洪水几乎完全消失的时候拍摄的。两幅图像的尺寸均为301301像素，灰度级为256，等效视数分别为10.89和9.26。图(c)为实际变化参考图，该参考图通过当地真实的陆地信息结合一些已知的图片信息而得到的，图中的白色像素表示两时相间发生变化的区域。l 实验结果图 （a) (b)图3.5 Bern数据结果图（a)差异图与(b)检测结果图图3.6 Bern 差异图直方图图3.7 Bern曲线拟合图（红色表示直方图曲线，蓝 色表示非变化类概率曲线） 采用第二章中用过的加拿大Ottawa地区两时相SAR图像为原始数据，用过本论文设计的算法的变化检测结果图。l 实验结果图 （a) (b)图3.8 Ottawa数据结果图（a)差异图，(b)检测结果图实验结果图/图3.9 Ottawa 差异图直方图图3.10 Ottawa曲线拟合图（红色表示直方图曲线，蓝色表示非变化类概率曲线）结果分析：Bern的检测结果图与原参考图相比，检测结果图漏检311个像素点，误检97个像素点。检测结果图接近参考结果图，但效果不是很理想，漏检数较多，误检数一般，需要改进。Ottawa检测结果图与原参考图相比，检测结果图漏检685个像素点，误检2317个像素点。检测结果图接近参考结果图，但效果不是很理想，漏检数较少，但错检数较多，需要改进。与参考图相比，本论文设计的方法虽然含有一定的噪声，但其对变化区域的检测更加准确。相对于一些整体拟合的策略方法，本方法能更好地拟合阈值产生区域的直方图，从而得到更佳的SAR图像变化检测结果图。但相对于第二章中提到的实际中常用的几种SAR图像变化检测方法，本方法的算法较为复杂，计算量大，其效果相对于图像差值法效果较好，但相对于对数比值法效果一般。通过对照一些文献中的设计的拟合方法，比如基于单边高斯对数模型的算法、基于log-normal分布模型的阈值算法等，本方法效果相对较差，可能在于此方法其中关键的两步，差异图的构造方法以