CN104361582B - 一种对象级高分辨率sar影像洪水灾害变化检测方法 - Google Patents

Abstract

本发明公开一种对象级高分辨率SAR影像洪水灾害变化检测方法，针对SAR影像中可能存在的与水体区域具有相似几何特征及波谱特征的虚假目标，造成大量的“伪变化”给洪水灾害变化检测造成的困难与干扰。因此，对每个时相影像进行轮廓波变换。在保持图像边缘特性的同时压制图像噪声，同时通过简单的样本训练选择最佳分解尺度，并在该尺度上通过分块直方图统计快速获取水体可能存在区域的标记点位置。进而采用基于标记点的分水岭分割以及区域合并策略获得水体的轮廓信息，并利用基于多种特征的判别规则进一步消除虚假目标的干扰。最后，根据配准结果直接比较多时相影像中所提取的水体轮廓，获得水体发生变化的区域。

Description

一种对象级高分辨率SAR影像洪水灾害变化检测方法

技术领域

本发明涉及一种对象级高分辨率SAR影像洪水灾害变化检测方法，属于高分辨SAR影像洪水灾害变化检测技术领域。

背景技术

洪水灾害是指由于降水导致的山洪爆发、堤坝溃坝、农田淹没、房屋损毁等现象，是目前世界上危害最为严重的自然灾害之一。遥感影像由于具有时效性强、覆盖范围广泛等优点，目前已成为了洪水灾害监测与评估的重要手段。由于水体对于014～215Lm范围内的电磁波吸收率明显高于其他地物，因此水体通常表现为面状连续分布且整体亮度较低的成片区域，这也是利用遥感影像进行洪水灾害变化检测的基础。

目前全球用于洪水灾害监测的卫星主要包括美国陆地资源卫星Landsat TM、法国SPOT系列资源卫星等光学遥感卫星，以及加拿大Radarsat SAR雷达卫星等。例如2012年巴拉圭利用空间分辨率为10米的SPOT 5卫星对洪水进行了监测，我国利用ERS－1/2和Radarsat影像监测了1998年发生在长江和嫩江流域的洪灾等。但是，由于传统资源卫星的轨道周期较长，因此难以获得动态的洪灾信息，时效性较差。同时，光学卫星对难以穿透云层的遮挡，在雨季很难获得可靠且清晰的可用影像。而合成孔径雷达SAR(SyntheticAperture Radar)具有全天候、全天时成像、能够穿透云层、几何分辨率与传感器的波长和高度无关等特性，从而能够有效揭示地面伪装及地貌结构，因此在洪水灾害监测中受到了高度重视。

近年来随着SAR影像空间分辨的不断提高(如分辨率为3m的意大利高分辨率雷达卫星COSMO-SkyMed SAR影像等)，为高分辨率SAR影像洪水灾害变化检测带来了新的机遇，同时也带来了新的挑战：空间分辨的提高使SAR影像中的地物纹理特征更为精细，层次更为丰富，但SAR影像固有的干斑噪声以及大量以纯象元形式出现的如小金属、小目标等噪声造成的干扰也更加明显，因此需要有效的滤波方法进行去噪，同时能够保持图像细节的特征；影像中与水体区域具有相似几何特征及波谱特征的虚假目标可能造成大量的“伪变化”；如何在背景复杂的场景中准确定位水体区域并准确提取水体的轮廓是变化检测成功的关键因素。

针对洪灾前后水体区域这一特定目标的变化检测，可首先对多时相影像单独进行分割，进而判断发生变化的区域，再通过分类提取出变化对象中属于水体的区域，从而获得最终检测结果。但这种方法在不会将水体区域作为一个独立的对象进行分析，因而在变化检测中忽略了水体区域独有的形状、光谱、面积等特征，难以获得可靠的水体区域变化检测结果。同时，变化检测产生的误差积累也会对分类精度造成影响。而在单一时相影像中先将水体提取作为独立的对象进行提取，再比较多时相影像间发生变化的水体区域，则能很好的避免这些问题和局限。由于变化检测的目的仅关注水体变化信息，因此分类问题可转化为如何准确提取影像中的水体区域这一关键问题。最后，只需根据配准结果直接比较多时相影像中提取的水体区域即可获得最终变化检测结果。因而准确提取高分辨率SAR影像中的水体轮廓、消除虚假目标造成的干扰是保证检测精度的关键因素。

针对高分辨率SAR影像水体提取及变化检测问题，学者们已开展了一些研究工作。目前主流的SAR影像水体检测方法主要分为三类:第一种是基于传统边缘检测的水体提取方法。如朱俊杰等采用块追踪算法初步确定水体边缘所在区域，进而采用活动轮廓(Snake)模型进一步精确定位水体边缘。但这种方法对噪声的鲁棒性较差，也容易受到虚假目标的干扰。第二种方法通常直接通过阈值法对图像进行二值化分割，从而获得水体及非水体区域。如戴光照等根据影像的灰度直方图采用阈值分割的方法实现了水体的提取。李金基等进一步提出首先对差值影像和对数比图像进行融合，进而采用一种基于T-分布模型改进的KI阈值法对融合影像进行分割，从而获得洪水灾害变化检测结果。这种方法原理简单，易于实现，但基于阈值的影像分割对高分辨率SAR影像中的细节信息不敏感，因而难以准确定位水体的边缘。第三种是基于支持向量机(SVM，Support Vector Machine)等分类器的检测方法。如程明跃等通过模糊加权支持向量机对样本图像的纹理进行训练，获得判别水体的决策函数，进而检测出SAR图像的水体区域，但这种基于分类器的检测算法复杂度较高，且样本获取的代价较大。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供一种对象级高分辨率SAR影像洪水灾害变化检测方法。针对SAR影像中可能存在的与水体区域具有相似几何特征及波谱特征的虚假目标可能造成大量的“伪变化”。因此，对每个时相影像进行轮廓波(Countourlet)变换。在保持图像边缘特性的同时压制图像噪声，同时通过简单的样本训练选择最佳分解尺度，并在该尺度上通过分块直方图统计快速获取水体可能存在区域的标记点位置。进而采用基于标记点的分水岭分割以及区域合并策略获得水体的轮廓信息，并利用基于多种特征的判别规则进一步消除虚假目标的干扰。最后，根据配准结果直接比较多时相影像中所提取的水体轮廓，获得水体发生变化的区域。

技术方案：一种对象级高分辨率SAR影像洪水灾害变化检测方法，主要由几何配准、基于轮廓波变换的噪声抑制及标记点提取、基于标记点的分水岭分割及区域合并、虚假目标消除及变化检测组成。

其中图像配准采用美国ITT视觉解决方案(Visual Information Solutions)的ENVI软件；

轮廓波变换主要包括两个步骤：

Step1：采用多尺度拉普拉斯金字塔分解影像，捕捉影像中的奇异点信息。

LP分解是多尺度分析方法，每次分解过程都首先采用低通滤波器对上一尺度进行滤波，进而进行下采样获得当前尺度的低频影像。对当前尺度低频影像做上采样操作，并进一步利用低通滤波进行平滑，最后与上一尺度低频影像相减可获得当前尺度下的高频影像。因此，LP分解没有对高频影像的方向进行划分，从而利于后续采用方向滤波器进一步提取高频分量的方向信息。

Step2：采用方向滤波器组(DFB，Directional Filter Bank)对高频影像进行方向滤波。方向滤波器能够将某一方向上不连续的点连接起来构成轮廓段，从而利用这种基结构来逼近原始影像；

针对SAR影像中固有的相干斑噪声，首先采用轮廓波变换对配准后的多时相SAR影像进行多尺度去噪。在轮廓波分解后的影像中，噪声对应较小的轮廓波系数，因此通常采用硬阈值方法对噪声加以区分。采用如下多尺度策略：

Step1：若原始影像尺寸为N×N像素，定义影像中像素的灰度值为f_i,j(i,j＝1,2,3....N)。定义分解层数M，进行方向为2^M的轮廓波变换。

Step2：利用公式(1)对分解后第m(m＝1,2...M)层的系数的方差δ²(m)进行估计，如下所示：

其中，C为常数。为第m层的轮廓波变换系数。

Step3：利用公式(2)计算第m层各高频子带的阈值，如下所示：

其中，L(m)为第m层分解系数的数量。

Step4：采用硬阈值方法对不同方向、不同子带的系数进行处理。

Step5：进行轮廓波反变换，获得去噪影像。

最佳轮廓波分解尺度选择：

通过计算所选取水体样本的局部方差均值，选择均值最小的尺度为最佳分解尺度，并认为在此尺度下水体内部的均质度最高，因此仅在此尺度下搜索可能存在水体的区域，并提取标记点。最佳轮廓波分解尺度选择过程如下：

在图像去噪中，进行M层轮廓波变换后，每层取K＝2^M个方向，则对图像分解得到M×K个高频子带的轮廓波系数W_m,k，k＝1,2,3,...K表示其对应的高频子带。分别利用公式(3)计算各尺度下各个波段的相位一致值PC_m,k：

其中，A_n为W_m,k傅里叶分解的n次谐波分量的幅度，φ_n(x)表示该相位偏移量分量在x处的局部相位，为所在点的傅里叶分量的加权平均相位。若所有傅里叶分量都具有一致的相位则该比值为1；反之该比值最小值为0。

对各尺度下各个波段的相位一致值PC_m,k，求梯度G_m,k。采用公式(4)融合个尺度下不同波段相位一致梯度值：

在原始图像中选取c个大小为d×d像素的水体区域作为样本，为了使个尺度下的样本具有相同的地物面积，在各个尺度所对应相位一致梯度图的相同位置用不同尺寸的窗口对水体样本进行采样，用公式(5)计算个尺度下水体的局部方差均值：

式中，σ² _m为尺度m下水体地物的局部方差均值；g×g像素为原始影像中样本区域在当前尺度下的对应区域的尺寸；表示m尺度下坐标为(x,y)的像元的相位一致梯度值；为采样样本的相位一致梯度均值。取σ_m ²值最小的尺度作为最佳分解尺度，并在此尺度下确定标记点位置。

基于分块直方图统计的标记点获取：

将待检测的低频图像进行分块处理，即将整幅影像等分为H×H像素大小的子图像，分别进行标记点提取。由于高分辨率SAR图像中水体的灰度值较低且波动范围小的特点，即反映为图像灰度的平均值和标准差较小，提出了基于分块直方图统计的标记点提取策略：

Step1：确定子图像尺寸。若窗口过大，窗口中可能包含种类众多的地物；如果窗口太小就会失去统计意义，同时会显著增加计算量。由于人工提取的样本区域能够充分反映水体所在区域的相关特征，因此选择子图像尺寸与所选样本在当前尺度下的对应区域尺寸一致，即令H＝g。定义任意一个子图像为I_i，且i＝1,2...U，U为子图像总数。

Step2：对I_i进行“灰度均值规则”判别：计算子图像的灰度均值μ_i，并与整幅影像的灰度均值μ进行比较。若满足μ_i≤μ，则将I_i标记为一个潜在的水体区域I′_i，进入下一步。否则不进行标记，重复Step2对下一个子图像进行判别。

Step3：对I′_i进行“相邻性规则”判别：若I′_i存在左侧水平相邻或者上方垂直相邻的子图像块，则以I′_i为中心在这两个方向上分别检测相邻的两个子图像，若有任一相邻子图像已被检测最终标记为水体区域块，则对当前I′_i不做任何标记，对下一块子图像进行判别。否则，将I′_i标记为I″_i，进入下一步。

Step4：对I″_i进行“直方图特征规则”判别：鉴于水体区域主要表现为亮度值较低的成片区域，计算I″_i的灰度直方图后，处于水体区域的子图像直方图应为左边单峰形状。满足如下条件的I″_i为水体区域：直方图峰值灰度级小于图像平均灰度级；峰值灰度级及其右边8个灰度级的像素数和占所在区域像素数比例大于60％。对满足条件的I″_i标记为I″′_i，取I″′_i的中心像素点作为标识水体主要区域的一个种子点；否则，不进行标记。重复Step2到Step4，遍历所有子图像，获得当前尺度下的标记点。

Step5：将当前尺度下提取的标记点映射到原始影像中的对应区域，取区域的中心作为最终的标记点位置。

基于标记点的分水岭分割及区域合并：

在利用标记点标记水体区域可能存在的位置后，对原始影像进行基于标记点的分水岭分割。首先计算原始影像的梯度影像，将梯度影像中的每个像素的像素值作为点的海拔高度值，将标记点作为积水盆地的底部。由于所提取的标记点仅标记了水体可能存在的区域，为避免欠分割现象，设定所提取标记点在梯度影像中的最大灰度值为阈值T_e，将所有灰度值小于T_e的像素同样作为一个盆地的底部。将周围的像素按照海拔从低到高(灰度值从小到大)逐个并入每个盆地底部所在的区域，构成相应的积水盆地。随着各个盆地的面积不断扩大，当两个盆地的水流相汇时的边界，即为分水岭，对应分割结果中相邻对象的边界。由于高分辨率SAR图像中地物种类众多，纹理特征丰富，在分割结果中可能存在过分割现象，尤其可能存在属于同一片水体的相邻区域，因此需要进一步的合并处理处理。

设定阈值T_d，由于参与合并的相邻区域应当具有相似的灰度分布，因此他们灰度均值向量的欧氏距离应位于阈值T_d的区间内。具体合并过程为：设分割结果中区域为R＝{R_i,i＝1,2,3,...,N}，N为区域总个数。各区域的灰度均值向量为μ＝(μ₁,μ₂,...,μ_N)，据此构建区域连接矩阵U＝{n_ij,i,j＝1,2,...,N}，若R_i和R_j相邻，则n_ij＝1，否则n_ij＝0。区域合并过程如下：

Step1：根据分割结果生成区域邻接图。

Step2：对每一个与R_i相邻的区域R_j，计算灰度均值向量间的欧式距离，若位于阈值T_d的区间内，则进行合并，并更新区域邻接图，重复Step2。否则，直接对下一个区域进行判别。

Step3:重复Step2，直到没有需要合并的区域为止。

由于所提出的标记点提取策略已经标记出所有水体可能存在的区域，因此经过区域合并处理后，分割结果仅保留包含标记点的部分区域，从而实现了水体区域的轮廓信息提取。

虚假目标消除：

通过以上步骤提取的水体区域依然可能存在一些虚假目标，为此进一步定义判别规则如下：

规则一：计算每个水体区域的平均灰度值，将大于整幅影像的灰度均值的区域判定为虚假目标。

规则二：生成水体区域的灰度直方图，利用上述中提取水体区域标记点时提出的“直方图特征规则”进行虚假目标判别。

规则三：计算每个水体区域的面积，设定阈值T_area，水体区域的面积应大于阈值T_area，否则为虚假目标。

符合以上规则的区域被判定为水体区域。最后，采用同样的策略分别提取多时相SAR影像中的水体区域，并根据配准结果对这些区域进行比较，获得最终洪水灾害变化检测结果。

有益效果：本发明采用基于标记点的分水岭分割及区域合并策略提取了水体区域的轮廓信息，并提出了虚假目标消除策略，实现了单一时相影像中水体轮廓的准确提取。最后，根据配准结果直接比较多时相影像中的水体区域，获得了洪水灾害变化检测结果。通过对两组单一时相高分辨率TerraSAR-X影像以及一组多时相高分辨率RADARSAT影像进行实验，可以得出如下结论：

1、所提出的基于最佳尺度的标记点获取策略是可行且有效的，能够准确标识出可能存在水体的区域，只有少数标记点受到了虚假目标的干扰，为进一步的水体提取奠定了良好基础。

2、所提出的基于分水岭的区域分割及合并以及虚假目标消除策略，能够准确提取水体区域的轮廓，并有效消除虚假目标的干扰。

3、以包含水体的区域为感兴趣对象，通过引入光谱、形状、面积、空间关系等多个特征进行标记点提取、影像分割、区域合并以及虚假目标消除，能够有效提取单一时相影像中水体区域及其轮廓细节特征。在多组实验中，本发明在目视分析及精度指标评价中都明显优于两种基于阈值分割的检测方法，可为高分辨率SAR影像洪水灾害变化检测提供一种可行且有效的解决方案。

附图说明

图1为本发明实施例的方法流程图；

图2为本发明实施例的拉普拉斯分解流程图；

图3为本发明实施例的轮廓波变换结构图；

图4中，(a)原始影像#3，(b)原始影像#4；

图5中，(a)影像#3水体区域，(b)影像#4水体区域；

图6为参考数据；

图7为本发明实施例方法的变化检测结果

图8为KI阈值法变化检测结果。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，本发明所提出方法主要由几何配准、基于轮廓波变换的噪声抑制及标记点提取、基于标记点的分水岭分割及区域合并、虚假目标消除及变化检测组成。

其中图像配准采用美国ITT视觉解决方案(Visual Information Solutions)的ENVI软件，不再赘述；

轮廓波采用轮廓段(Contour Segment)的基结构来逐渐逼近对象的轮廓，由于定义了能够根据尺度变化而改变长宽比的长方形结构作为基的支撑区间，因而具有各向异性及方向性。轮廓波变换主要包括两个步骤：

Step1：采用多尺度拉普拉斯金字塔分解影像，捕捉影像中的奇异点信息。

LP分解是多尺度分析方法，每次分解过程都首先采用低通滤波器对上一尺度进行滤波，进而进行下采样获得当前尺度的低频影像。对当前尺度低频影像做上采样操作，并进一步利用低通滤波进行平滑，最后与上一尺度低频影像相减可获得当前尺度下的高频影像。因此，LP分解没有对高频影像的方向进行划分，从而利于后续采用方向滤波器进一步提取高频分量的方向信息。分解实现流程如图2所示。

其中，l为当前尺度的低频影像，l+1为分解后获得的低频影像，k为输出的带通影像，H为分解滤波器，G为综合滤波器。

Step2：采用方向滤波器组(DFB，Directional Filter Bank)对高频影像进行方向滤波。方向滤波器能够将某一方向上不连续的点连接起来构成轮廓段，从而利用这种基结构来逼近原始影像；轮廓波变换的结构示意图如图3所示。

针对SAR影像中固有的相干斑噪声，首先采用轮廓波变换对配准后的多时相SAR影像进行多尺度去噪。在轮廓波分解后的影像中，噪声对应较小的轮廓波系数，因此通常采用硬阈值方法对噪声加以区分。采用如下多尺度策略：

Step1：若原始影像尺寸为N×N像素，定义影像中像素的灰度值为f_i,j(i,j＝1,2,3....N)。定义分解层数M，进行方向为2^M的轮廓波变换。

Step2：利用公式(1)对分解后第m(m＝1,2...M)层的系数的方差δ²(m)进行估计，如下所示：

其中，C为常数。为第m层的轮廓波变换系数。

Step3：利用公式(2)计算第m层各高频子带的阈值，如下所示：

其中，L(m)为第m层分解系数的数量。

Step4：采用硬阈值方法对不同方向、不同子带的系数进行处理。

Step5：进行轮廓波反变换，获得去噪影像。

最佳轮廓波分解尺度选择：

在利用轮廓波进行图像去噪的同时，进一步利用轮廓波分解后的低频影像对水体区域可能存在的位置进行标记。在轮廓波分解后的不同尺度低频影像中，地物特征的详细程度会发生变化。低频影像尺度越大，图像尺寸越小，包含的细节信息也逐渐减少。同时，一些小目标、空洞及孤立点等可能会消失，而如水体、道路等轮廓特征明显、面积较大的典型地物却更加突出。因此，选择合适的低频影像标记水体可能存在的区域，不但能够减少运算量，同时利于提高所获取标记点的准确性。考虑到水体区域在高分辨率SAR图像中通常表现为灰度较为稳定的成片区域，利用局部方差对不同尺度下水体样本的光谱特征进行统计将具有很好的可比性。本发明通过计算所选取水体样本的局部方差均值，选择均值最小的尺度为最佳分解尺度，并认为在此尺度下水体内部的均质度最高，因此仅在此尺度下搜索可能存在水体的区域，并提取标记点。最佳轮廓波分解尺度选择过程如下：

在图像去噪中，进行M层轮廓波变换后，每层取K＝2^M个方向，则对图像分解得到M×K个高频子带的轮廓波系数W_m,k，k＝1,2,3,...K表示其对应的高频子带。分别利用公式(3)计算各尺度下各个波段的相位一致值PC_m,k：

其中，A_n为W_m,k傅里叶分解的n次谐波分量的幅度，φ_n(x)表示该相位偏移量分量在x处的局部相位，为所在点的傅里叶分量的加权平均相位。若所有傅里叶分量都具有一致的相位则该比值为1；反之该比值最小值为0。

对各尺度下各个波段的相位一致值PC_m,k，求梯度G_m,k。采用公式(4)融合个尺度下不同波段相位一致梯度值：

在原始图像中选取c个大小为d×d像素的水体区域作为样本，为了使个尺度下的样本具有相同的地物面积，在各个尺度所对应相位一致梯度图的相同位置用不同尺寸的窗口对水体样本进行采样，用公式(5)计算个尺度下水体的局部方差均值：

式中，σ² _m为尺度m下水体地物的局部方差均值；g×g像素为原始影像中样本区域在当前尺度下的对应区域的尺寸；表示m尺度下坐标为(x,y)的像元的相位一致梯度值；为采样样本的相位一致梯度均值。取σ_m ²值最小的尺度作为最佳分解尺度，并在此尺度下确定标记点位置。

基于分块直方图统计的标记点获取：

经过轮廓波变换后的低频子图像保留了原始图像的大部分稳定信息，且分辨率相对于原始图像小，因此即使在最佳尺度对应的低频子图像中进行全局搜索，搜索次数也会有效减少。为了进一步减少计算量本发明将待检测的低频图像进行分块处理，即将整幅影像等分为H×H像素大小的子图像，分别进行标记点提取。由于高分辨率SAR图像中水体的灰度值较低且波动范围小的特点，即反映为图像灰度的平均值和标准差较小，提出了基于分块直方图统计的标记点提取策略：

Step1：确定子图像尺寸。若窗口过大，窗口中可能包含种类众多的地物；如果窗口太小就会失去统计意义，同时会显著增加计算量。由于人工提取的样本区域能够充分反映水体所在区域的相关特征，因此选择子图像尺寸与所选样本在当前尺度下的对应区域尺寸一致，即令H＝g。定义任意一个子图像为I_i，且i＝1,2...U，U为子图像总数。

Step2：对I_i进行“灰度均值规则”判别：计算子图像的灰度均值μ_i，并与整幅影像的灰度均值μ进行比较。若满足μ_i≤μ，则将I_i标记为一个潜在的水体区域I′_i，进入下一步。否则不进行标记，重复Step2对下一个子图像进行判别。

Step3：对I′_i进行“相邻性规则”判别：若I′_i存在左侧水平相邻或者上方垂直相邻的子图像块，则以I′_i为中心在这两个方向上分别检测相邻的两个子图像，若有任一相邻子图像已被检测最终标记为水体区域块，则对当前I′_i不做任何标记，对下一块子图像进行判别。否则，将I′_i标记为I″_i，进入下一步。

Step4：对I″_i进行“直方图特征规则”判别：鉴于水体区域主要表现为亮度值较低的成片区域，计算I″_i的灰度直方图后，处于水体区域的子图像直方图应为左边单峰形状。满足如下条件的I″_i为水体区域：直方图峰值灰度级小于图像平均灰度级；峰值灰度级及其右边8个灰度级的像素数和占所在区域像素数比例大于60％。对满足条件的I″_i标记为I″′_i，取I″′_i的中心像素点作为标识水体主要区域的一个种子点；否则，不进行标记。重复Step2到Step4，遍历所有子图像，获得当前尺度下的标记点。

Step5：将当前尺度下提取的标记点映射到原始影像中的对应区域，取区域的中心作为最终的标记点位置。

基于标记点的分水岭分割及区域合并：

在利用标记点标记水体区域可能存在的位置后，对原始影像进行基于标记点的分水岭分割。首先计算原始影像的梯度影像，将梯度影像中的每个像素的像素值作为点的海拔高度值，将标记点作为积水盆地的底部。由于所提取的标记点仅标记了水体可能存在的区域，为避免欠分割现象，设定所提取标记点在梯度影像中的最大灰度值为阈值T_e，将所有灰度值小于T_e的像素同样作为一个盆地的底部。将周围的像素按照海拔从低到高(灰度值从小到大)逐个并入每个盆地底部所在的区域，构成相应的积水盆地。随着各个盆地的面积不断扩大，当两个盆地的水流相汇时的边界，即为分水岭，对应分割结果中相邻对象的边界。由于高分辨率SAR图像中地物种类众多，纹理特征丰富，在分割结果中可能存在过分割现象，尤其可能存在属于同一片水体的相邻区域，因此需要进一步的合并处理处理。

设定阈值T_d，由于参与合并的相邻区域应当具有相似的灰度分布，因此他们灰度均值向量的欧氏距离应位于阈值T_d的区间内。具体合并过程为：设分割结果中区域为R＝{R_i,i＝1,2,3,...,N}，N为区域总个数。各区域的灰度均值向量为μ＝(μ₁,μ₂,...,μ_N)，据此构建区域连接矩阵U＝{n_ij,i,j＝1,2,...,N}，若R_i和R_j相邻，则n_ij＝1，否则n_ij＝0。区域合并过程如下：

Step1：根据分割结果生成区域邻接图。

Step2：对每一个与R_i相邻的区域R_j，计算灰度均值向量间的欧式距离，若位于阈值T_d的区间内，则进行合并，并更新区域邻接图，重复Step2。否则，直接对下一个区域进行判别。

Step3:重复Step2，直到没有需要合并的区域为止。

由于所提出的标记点提取策略已经标记出所有水体可能存在的区域，因此经过区域合并处理后，分割结果仅保留包含标记点的部分区域，从而实现了水体区域的轮廓信息提取。

虚假目标消除：

通过以上步骤提取的水体区域依然可能存在一些虚假目标，为此进一步定义判别规则如下：

规则一：计算每个水体区域的平均灰度值，将大于整幅影像的灰度均值的区域判定为虚假目标。

规则二：生成水体区域的灰度直方图，利用上述中提取水体区域标记点时提出的“直方图特征规则”进行虚假目标判别。

规则三：计算每个水体区域的面积，设定阈值T_area，水体区域的面积应大于阈值T_area，否则为虚假目标。

符合以上规则的区域被判定为水体区域。最后，采用同样的策略分别提取多时相SAR影像中的水体区域，并根据配准结果对这些区域进行比较，获得最终洪水灾害变化检测结果。

实验结果与分析：

为了进一步验证所提出的算法在洪水灾害变化检测中的性能，选择了一组双时相加拿大渥太华(Ottawa)地区1997年洪灾前后的高分辨率RADARSAT影像进行实验，并与李金基等人提出的洪水灾害变化检测算法(下文中简称为“KI阈值法”)进行了比较和分析。KI阈值法结合了插值法与对数比的优点，实现了差分影像的增强，并在此基础上利用一种基于T-分布模型改进的KI阈值法实现了对差分影像的阈值分割获得了洪水灾害变化结果，取得了良好的效果。

水体区域变化检测实验结果与分析：

为进一步验证本文方法的在洪水灾害变化检测中的效果，第三组数据采用两幅加拿大渥太华(Ottawa)地区的洪水灾害前后的高分辨率RADARSAT影像进行实验，空间分辨率为8.9m，影像尺寸为256×256像素，并与KI阈值法进行了比较实验。进行几何配准后的影像#3，影像#4如图4(a)，图4(b)所示。其中影像#3采集时间为1997年5月，影像#4采集时间为1997年8月，两幅影像反映的为1997年5月夏季暴雨导致的洪水灾害前后水体区域变化。

通过上图可以看出，场景中存在大面积的水体、湖泊以及零星分布的小片水体区域。设定T_area＝200像素，T_d＝[-2,2]，其他参数设定均与前两组实验相同。在KI阈值法中，设定Lee滤波参数N_L＝3，融合参数为N_F＝7，S_T＝0.7。

采用本发明方法提取的水体轮廓如图5(a)，图5(b)所示，其中白色区域为水体。

为评价检测精度，通过目视解译提取了水体发生变化的区域作为参考数据，如图6所示。基于配准结果比较影像#3和影像#4中的水体区域，获得水体发生变化的区域如图7所示。图8为KI阈值法提取的水体变化区域。

通过目视分析可以看出，两种算法都提取出了洪水灾害前后水体变化的主要区域，而本文提出算法的检测结果最为接近图6的参考影像。KI阈值法尽管采用了基于差分影像与对比影像的图像增强预处理，但基于阈值分割的变化检测策略仍然难以避免图像细节特征的损失，从而无法提取水体边缘细微的变化。进一步对两种算法的实验结果进行精度评价如表1所示：

表1精度评价

通过表1可以看出，KI阈值法在预处理过程中增加了图像的对比度，因而在抑制误检率方面取得了良好的效果。但本发明方法在总体精度等评价指标中明显优于KI阈值法，与目视分析结果一致。

Claims (4)

1.一种对象级高分辨率SAR影像洪水灾害变化检测方法，其特征在于，主要由几何配准、基于轮廓波变换的噪声抑制及标记点提取、基于标记点的分水岭分割及区域合并、虚假目标消除及变化检测组成；

其中图像配准采用美国ITT视觉解决方案的ENVI软件；

轮廓波变换主要包括两个步骤：

Step1：采用多尺度拉普拉斯金字塔分解影像，捕捉影像中的奇异点信息；

每次分解过程都首先采用低通滤波器对上一尺度进行滤波，进而进行下采样获得当前尺度的低频影像；对当前尺度低频影像做上采样操作，并进一步利用低通滤波进行平滑，最后与上一尺度低频影像相减可获得当前尺度下的高频影像；

Step2：采用方向滤波器组对高频影像进行方向滤波；

针对SAR影像中固有的相干斑噪声，首先采用轮廓波变换对配准后的多时相SAR影像进行多尺度去噪；在轮廓波分解后的影像中，噪声对应较小的轮廓波系数，采用硬阈值方法对噪声加以区分；采用如下多尺度策略：

Step1：若原始影像尺寸为N×N像素，定义影像中像素的灰度值为f_i,j，i,j＝1,2,3....N；定义分解层数M，进行方向为2^M的轮廓波变换；

Step2：利用公式(1)对分解后第m层的系数的方差δ²(m)进行估计，其中m＝1,2...M，如下所示：

${{\mathrm{\δ}}_{f_{i,j}}}^2\left(m\right)=\frac{Med\left(\right|W_{f_{i,j}}\left(m\right)\left|\right)}{C}---\left(1\right)$

其中，C为常数；为第m层的轮廓波变换系数；

Step3：利用公式(2)计算第m层各高频子带的阈值，如下所示：

$T\left(m\right)={{\mathrm{\δ}}_{f_{i,j}}}^2\left(m\right)\sqrt{2lgL\left(m\right)}---\left(2\right)$

其中，L(m)为第m层分解系数的数量；

Step4：采用硬阈值方法对不同方向、不同子带的系数进行处理；

Step5：进行轮廓波反变换，获得去噪影像；

最佳轮廓波分解尺度选择：

通过计算所选取水体样本的局部方差均值，选择均值最小的尺度为最佳分解尺度，并认为在此尺度下水体内部的均质度最高，因此仅在此尺度下搜索可能存在水体的区域，并提取标记点；最佳轮廓波分解尺度选择过程如下：

在图像去噪中，进行M层轮廓波变换后，每层取K＝2^M个方向，则对图像分解得到M×K个高频子带的轮廓波系数W_m,k，k＝1,2,3,...K表示其对应的高频子带；分别利用公式(3)计算各尺度下各个波段的相位一致值PC_m,k：

${\mathrm{PC}}_{m,k}=\frac{{\Σ}_n{A}_n\left(x\right)cos\left({\mathrm{\φ}}_n\right(x)-\stackrel{\‾}{\mathrm{\φ}}(x\left)\right)}{\underset{n}{\Σ}{A}_n\left(x\right)}---\left(3\right)$

其中，A_n为W_m,k傅里叶分解的n次谐波分量的幅度，φ_n(x)表示该相位偏移量分量在x处的局部相位，为所在点的傅里叶分量的加权平均相位；若所有傅里叶分量都具有一致的相位则公式(3)的比值为1；反之公式(3)的比值最小值为0；

对各尺度下各个波段的相位一致值PC_m,k，求梯度G_m,k；采用公式(4)融合个尺度下不同波段相位一致梯度值：

${\tilde{G}}_m=max(G_{m1},G_{m2},...,G_{mK})---\left(4\right)$

在原始图像中选取c个大小为d×d像素的水体区域作为样本，为了使个尺度下的样本具有相同的地物面积，在各个尺度所对应相位一致梯度图的相同位置用不同尺寸的窗口对水体样本进行采样，用公式(5)计算个尺度下水体的局部方差均值：

${{\mathrm{\σ}}^2}_m=\frac{1}{c}\Σ\frac{1}{g\×g}\underset{x=1}{\overset{g}{\Σ}}\underset{y=1}{\overset{g}{\Σ}}{\[{\tilde{G}}_m(x,y)-{\stackrel{\‾}{\tilde{G}}}_m\]}^2---\left(5\right)$

式中，σ² _m为尺度m下水体地物的局部方差均值；g×g像素为原始影像中样本区域在当前尺度下的对应区域的尺寸；表示m尺度下坐标为(x,y)的像元的相位一致梯度值；为采样样本的相位一致梯度均值；取σ_m ²值最小的尺度作为最佳分解尺度，并在此尺度下确定标记点位置。

2.如权利要求1所述的对象级高分辨率SAR影像洪水灾害变化检测方法，其特征在于，将待检测的低频图像进行分块处理，即将整幅影像等分为H×H像素大小的子图像，分别进行标记点提取；由于高分辨率SAR图像中水体的灰度值较低且波动范围小的特点，即反映为图像灰度的平均值和标准差较小，采用基于分块直方图统计的标记点提取策略：

Step1：确定子图像尺寸；选择子图像尺寸与所选样本在当前尺度下的对应区域尺寸一致，即令H＝g，定义任意一个子图像为I_i，且i＝1,2...U，U为子图像总数；

Step2：对I_i进行“灰度均值规则”判别：计算子图像的灰度均值μ_i，并与整幅影像的灰度均值μ进行比较；若满足μ_i≤μ，则将I_i标记为一个潜在的水体区域I_i′，进入下一步；否则不进行标记，重复Step2对下一个子图像进行判别；

Step3：对I_i′进行“相邻性规则”判别：若I_i′存在左侧水平相邻或者上方垂直相邻的子图像块，则以I_i′为中心在这两个方向上分别检测相邻的两个子图像；若有任一相邻子图像已被检测最终标记为水体区域块，则对当前I_i′不做任何标记，对下一块子图像进行判别；否则，将I_i′标记为I_i″，进入下一步；

Step4：对I_i″进行“直方图特征规则”判别：鉴于水体区域主要表现为亮度值较低的成片区域，计算I_i″的灰度直方图后，处于水体区域的子图像直方图应为左边单峰形状；满足如下条件的I_i″为水体区域：直方图峰值灰度级小于图像平均灰度级；峰值灰度级及其右边8个灰度级的像素数和占所在区域像素数比例大于60％；对满足条件的I_i″标记为I_i″′，取I_i″′的中心像素点作为标识水体主要区域的一个种子点；否则，不进行标记；重复Step2到Step4，遍历所有子图像，获得当前尺度下的标记点；

Step5：将当前尺度下提取的标记点映射到原始影像中的对应区域，取区域的中心作为最终的标记点位置。

3.如权利要求1所述的对象级高分辨率SAR影像洪水灾害变化检测方法，其特征在于，基于标记点的分水岭分割及区域合并：

在利用标记点标记水体区域可能存在的位置后，对原始影像进行基于标记点的分水岭分割；首先计算原始影像的梯度影像，将梯度影像中的每个像素的像素值作为点的海拔高度值，将标记点作为积水盆地的底部；设定所提取标记点在梯度影像中的最大灰度值为阈值T_e，将所有灰度值小于T_e的像素同样作为一个盆地的底部；将周围的像素按照海拔从低到高逐个并入每个盆地底部所在的区域，构成相应的积水盆地；随着各个盆地的面积不断扩大，当两个盆地的水流相汇时的边界，即为分水岭，对应分割结果中相邻对象的边界；

设定阈值T_d，具体合并过程为：设分割结果中区域为R＝{R_i,i＝1,2,3,...,N}，N为区域总个数；各区域的灰度均值向量为μ＝(μ₁,μ₂,...,μ_N)，据此构建区域连接矩阵U＝{n_ij,i,j＝1,2,...,N}，若R_i和R_j相邻，则n_ij＝1，否则n_ij＝0；区域合并过程如下：

Step1：根据分割结果生成区域邻接图；

Step2：对每一个与R_i相邻的区域R_j，计算灰度均值向量间的欧式距离，若位于阈值T_d的区间内，则进行合并，并更新区域邻接图，重复Step2；否则，直接对下一个区域进行判别；

Step3:重复Step2，直到没有需要合并的区域为止；

由于所提出的标记点提取策略已经标记出所有水体可能存在的区域，因此经过区域合并处理后，分割结果仅保留包含标记点的部分区域，从而实现了水体区域的轮廓信息提取。

4.如权利要求3所述的对象级高分辨率SAR影像洪水灾害变化检测方法，其特征在于，提取的水体区域依然可能存在一些虚假目标，为此进一步定义判别规则如下：

规则一：计算每个水体区域的平均灰度值，将大于整幅影像的灰度均值的区域判定为虚假目标；

规则二：生成水体区域的灰度直方图，利用“直方图特征规则”进行虚假目标判别；

规则三：计算每个水体区域的面积，设定阈值T_area，水体区域的面积应大于阈值T_area，否则为虚假目标；

符合以上规则的区域被判定为水体区域；最后，采用同样的策略分别提取多时相SAR影像中的水体区域，并根据配准结果对这些区域进行比较，获得最终洪水灾害变化检测结果。