基于灰度统计及区域编码的SAR图像水域识别方法

摘要

本发明公开了一种基于灰度统计及区域编码的SAR图像水域识别方法，它属于SAR图像处理技术领域。主要解决了现有技术难以区分水域与低灰度的植被和路网的问题。其实现步骤为：先对图像进行像素灰度归一化处理，然后利用本发明确定的灰度经验阈值将图像像素归类为目标区和非水域，再计算归类处理后图像中的连通区域面积以消除小面积的非水域，最后通过跑长读取剩余目标区的编码，根据读取的区域编码消除间断性的植被和路网非水域，最终识别出真实水域。本发明能够有效地解决SAR图像水域识别对噪声敏感，低灰度植被和路网难以区分的问题，具有良好快速的识别效果，可用于河流和湖泊的识别及水上目标的检测。

权利要求

1.一种基于灰度统计及区域编码的SAR图像水域识别方法，包括如下步骤：(1)对输入图像采用像素灰度归一化进行处理；

(2)依据灰度经验阈值，将归一化后的图像像素归类为目标区和非水域：(2a)对于像素灰度分布在0～1之间的情况，取经验阈值为0.255；对于像素灰度分布在0～255之间的情况，取经验阈值为65；

(2b)将像素大于经验阈值的判为目标区其像素值置1，小于阈值的判为非水域其像素值置0；

(3)计算目标区中连通区域面积，消除目标区内的非水域：(3a)统计目标区内连通区域像素点的个数为每个区域对应的面积；

(3b)用区域面积总和除以区域总个数得到一个总区域面积均值；

(3c)将目标区内每个区域的面积小于面积均值的区域作为非水域，并将其像素值置0予以消除；

(4)在每个区域的最长方向的1/2和1/3处，通过跑长读取目标区中每个区域的编码方式，利用区域编码的间断性或连续性，对目标区中剩余的连通区域进行筛选，消除间断性植被和路网的非水域，得到最终的水域识别结果，所述的通过跑长读取目标区中每个区域的编码方式，按如下步骤进行：首先，将目标区中每个区域的最小坐标点和最大坐标点的横、纵坐标之差的较大者作为该区域的长轴；

其次，在垂直于该区域长轴的1/3和2/3处，按照0和1出现的先后次序分别记录每次

0和1出现的次数及每次出现的长度；

所述的利用区域编码的间断性或连续性，对目标区中剩余的区域进行筛选，是在得到的区域编码中，将0和1交替出现次数大于9的区域，或者1的出现长度均小于10的区域判为间断性非水域；将0和1交替出现次数小于9的区域，或者1的出现次数大于25的区域判为连续性水域。

说明书

基于灰度统计及区域编码的SAR图像水域识别方法

技术领域

[0001] 本发明属于图像处理技术领域，具体涉及合成孔径雷达(Synthetic Aperture Radar)SAR图像水域识别，本发明用于水域的分割和识别，洪涝灾害检测及多时相、多传感器间同一地物的图像配准。

背景技术

[0002] 水域信息的提取是地理信息测绘更新、水利资源调查、洪涝灾害的监测、水上目标如桥梁、船舶等识别、城市湿地保护、多传感器或多时相遥感图像配准等工作必不可少的基础步骤。作为一种主动的微波遥感雷达，SAR所具有的全天时、全天候、穿透性强等优点，使得SAR图像成为地面目标监测的研究重点和热点。

[0003] SAR图像的水域识别研究中已有多种方法，如阈值法、聚类方法、区域增长、主动轮廓模型(Active contour Model，也称蛇算法Snake Model)等。通过区分不同目标的灰度阈值可以划分出图像中的目标和背景。此方法计算简单，运算效率高，但由于难以区分与水域灰度相近的非水域如植被、阴影、路网等，在水域识别中较少单独使用。最小经验熵法是阈值法的变形，依据光滑、均匀的区域熵很低。但是该方法对低灰度的阴影和植被不能有效区分。聚类方法和区域增长方法通过选择聚类中心或种子采用一定的隶属规则迭代更新直到停止条件。此类方法中初始点的选取对最终结果影响较大，并且运算量大，计算时间长，同样对水域区中的植被、阴影、路网等非同质区的效果较差。蛇算法通常先对SAR图像滤波，然后用蛇算法勾画出真实的水域轮廓的方法。然而蛇算法在迭代计算过程中对于噪声非常敏感，常常会收敛到部分噪声点，甚至出现不能收敛的现象。对于具有乘性噪声的SAR图像，蛇算法的水域识别效果更多地依赖于其前一步的滤波结果，而不同形式的滤波对SAR图像的水域真实边缘轮廓都会产生改变，最终识别结果并不理想；另外蛇算法的迭代运算复杂、计算量大，运算速度非常慢。另有研究者提出一种蛇算法的变形方法，即将水域轮廓上每点的坐标及方向表示为一组离散的点集序列，由贝叶斯准则判断该点周围的像素点，使能量最小并不断更新点集，最终得到水域轮廓。该方法中初始点集可以人工或自动选取，但是初始点集的位置影响最终的分割结果，与蛇算法存在同样的问题，方法的实用性差。

[0004] 因此人们提出实现水域识别的最佳数据源是时间和空间分辨率都较高的遥感数据，如SPOT、Landsat等，或者是SAR图像和光学遥感影像结合的方法。但是实际中要求同时具备同一地区的SAR图像、光学遥感影像和地理信息是有一定的困难；其次SPOT、Landsat等光学遥感影像的价格昂贵，用户难以长期忍受；再者，光学遥感影像易受云层、雨雾等影响，对天气条件的要求较高，影响水域识别的效果。

发明内容

[0005] 本发明的目的在于针对上述现有技术在SAR图像水域识别中，存在乘性噪声影响的植被和路网难区分问题，提出了一种基于灰度统计及区域编码的SAR图像水域识别方法，提高水域识别的精度。

[0006] 为实现上述目的，本发明的技术方案包括如下步骤：

[0007] (1)对输入图像采用像素灰度归一化进行处理；

[0008] (2)依据灰度经验阈值，将归一化后的图像像素归类为目标区和非水域：

[0009] (2a)对于像素灰度分布在0～1之间的情况，取经验阈值为0.255；对于像素灰度分布在0～255之间的情况，取经验阈值为65；

[0010] (2b)将像素大于经验阈值的判为目标区其像素值置1，小于阈值的判为非水域其像素值置0；

[0011] (3)计算目标区中连通区域面积，消除目标区内的非水域：

[0012] (3a)统计目标区内连通区域像素点的个数为每个区域对应的面积；

[0013] (3b)用区域面积总和除以区域总个数得到一个总区域面积均值；

[0014] (3c)将目标区内每个区域的面积小于面积均值的区域作为非水域，并将其像素值置0予以消除；

[0015] (4)在每个区域的最长方向的1/2和1/3处，通过跑长读取目标区中每个区域的编码方式，利用区域编码的间断性或连续性，对目标区中剩余的连通区域进行筛选，消除间断性植被和路网的非水域，得到最终的水域识别结果，

[0016] 所述的通过跑长读取目标区中每个区域的编码方式，按如下步骤进行：

[0017] 首先，将目标区中每个区域的最小坐标点和最大坐标点的横、纵坐标之差的较大者作为该区域的长轴；

[0018] 其次，在垂直于该区域长轴的1/3和2/3处，按照0和1出现的先后次序分别记录每次0和1出现的次数及每次出现的长度；

[0019] 所述的利用区域编码的间断性或连续性，对目标区中剩余的区域进行筛选，是在得到的区域编码中，将0和1交替出现次数大于9的区域，或者1的出现长度均小于10的区域判为间断性非水域；将0和1交替出现次数小于9的区域，或者1的出现次数大于25的区域判为连续性水域。

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

[0021] A、本发明由于采用跑长读取目标区中每个区域的编码方式，因而不受SAR图像乘性噪声的影响，克服了去噪对SAR图像的影响。

[0022] B、本发明由于采用少量的加减运算，易于硬件实现，识别效果良好。

[0023] 仿真试验证明，本发明对SAR图像和满足类似微波水域反射特性的光学遥感影像均有效，其Kappa系数达到0.91，平均运行速度约3.53s。

附图说明

[0024] 图1是本发明实现流程图；

[0025] 图2是本发明读取跑长编码过程示例图，其中：

[0026] 图2(a)是原始SAR图像，

[0027] 图2(b)是本发明的阈值分割结果图，

[0028] 图2(c)是本发明读取跑长矩阵的示例图；

[0029] 图3是本发明对SAR图像水域识别的结果图，其中：

[0030] 图3(a)是原始SAR图像，

[0031] 图3(b)是本发明阈值分割消除小面积区域后的结果图，

[0032] 图3(c)是本发明的最终水域识别结果图；

[0033] 图4是本发明对光学遥感影像水域识别的结果图，其中：

[0034] 图4(a)是原光学遥感影像，

[0035] 图4(b)是本发明阈值分割消除小面积区域后的结果图，

[0036] 图4(c)是本发明的最终水域识别结果图。

具体实施方式

[0037] 参照图1，本发明的具体实现步骤如下：

[0038] 步骤1，对输入图像进行归一化处理。

[0039] 对输入图像进行灰度归一化，将图像的灰度从实际灰度分布数规范化到0～1或者0～255之间，使图像灰度值分布一致。经此处理后的SAR图像水域低灰度特点保留，同时水域和非水域分界更显著。

[0040] 步骤2，依据经验分割阈值，将归一化后的图像划分为目标区和非水域。

[0041] 由于微波对光滑水域产生镜面反射，使得回波信号强度非常弱，所成SAR图像的相应连通区域的像素灰度值在较低的范围内分布；而微波对其他粗糙地物产生反射和散射，使得回波信号强度强弱不均且总体较强，使所成SAR图像的相应区域像素灰度值总体在较高的范围内分布。经过大量SAR图像水域区的灰度值统计分析，发现水域像素灰度分布范围为图像整体灰度分布范围的最低值到低25％处。若像素灰度分布在0～1之间，则水域像素灰度分布在0～0.25之间；若像素灰度分布在0～255之间，则水域像素灰度分布在0～64之间。据此确定水域灰度的经验阈值范围为整体灰度分布的低21％～26％范围内。在进行仿真实验时对于像素灰度分布在0～255之间的情况，一般取为65。

[0042] 本发明对归一化后的图像像素大于经验阈值的判为目标区，像素值置1；小于阈值的判为非水域，像素值置0。

[0043] 步骤3，计算目标区连通区域的面积，消除小面积非水域。

[0044] 阈值分割后的目标区内包含多个独立的连通区域，统计目标区每个独立的连通区域所包含的像素数作为对应独立区域的面积。将总面积除以总的区域个数，得到面积均值。将面积小于面积均值的区域判为非水域，并将其像素值置0，以消除小面积非水域。

[0045] 步骤4，跑长读取区域编码，消除间断性植被和路网。

[0046] 消除小面积非水域后的目标区中仍然包括植被和路网，本发明通过跑长读取每个区域的编码来消除植被和路网。粗糙地物的SAR回波信号强度强弱不均且总体较强，使所成SAR图像的相应区域像素灰度值总体在21％～26％以上范围内分布，但植被和路网的像素灰度值有一些分布在水域像素范围21％～26％以下范围内，有些分布在水域像素范围21％～26％以上范围内，造成植被和路网区域呈现出孔洞状，如图2(b)中的B区所示的植被区呈现的空洞状。相反，在阈值分割及消除小面积区域后，水域区呈现大面积的连通一致性，如图2(b)中的A区所示。

[0047] 水域区跑长读取的编码为连续不间断的1，植被和路网跑长读取的编码为0和1交替间断出现如图2(b)中以白色表示1，黑色表示0。例如，对于图2(b)中B区的非水域连通区读取其跑长编码，读取位置为图2(c)所示，左边一条白色线读取的编码为：

[0048] [0 0 1 1 1 0 1 1 1 1 0 0 0 0 0 0 0 0 1 1......0 0 1 1 0 0 1 1 1 0]，[0049] 图2(b)中A区域读取的编码为：

[0050] [0 1 1 1 1 1 1......1 1 1 1 1 1 1 0 0 0 0]，

[0051] 图2(b)中B和图2(b)中C两个不同区域跑长读取的编码不同，即

[0052] 图2(b)中的B区是0和1的频繁交替变化，而图2(b)中的C区是连续的1。

[0053] 跑长读取编码的具体步骤如下：

[0054] 1.将经过阈值分割和消除小面积区后的目标区像素值置1，非水域像素置0。

[0055] 2.用数字对目标区中每个独立的连通区域编号1：NR，连通区域内的像素点编号数字相同。

[0056] 3.从编号1开始查找其对应区域内的最小坐标点和最大坐标点之差，将差大者作为该区域的长方向。在垂直于长方向的1/3和2/3处，确定读取跑长编码的位置。其中分别再1/3和2/3处分别读取是通过实验得出的，仅仅在1/3并不能完全通过跑长矩阵体现对应区域的特点，而增加更多的读取跑长编码位置则大大浪费了程序运行时间。

[0057] 4.在跑长位置统计每次0和1交替出现的次数，0和1每次出现的长度，将0和1交替出现次数大于9的区域，或者1的出现长度均小于10的区域判为间断性非水域，将其对应像素值置0予以消除；将0和1交替出现次数小于9的区域，或者1的出现次数大于25的区域判为连续性水域对应像素值保持为1。

[0058] 5.重复步骤2至步骤4，直到遍历到编号NR的所有区域结束。

[0059] 本发明的效果可以通过仿真实验具体说明：

[0060] 1.实验条件

[0061] 实验所用微机CPU为Intel Pentiun43.0GHz内存1GB，实验平台为Matlab 7.0.1。实验图像数据为包含水域的SAR图像及满足微波水域反射特性的光学遥感影像。

[0062] 实验数据为搜集的35幅包含水域的遥感图像，分别以1～35编号，如表1所示。其中第1到第11幅和第25幅图像为Corona卫星KH-4拍摄的光学遥感图像；第19、20和
21幅为Landsat 5TM影像；第26和35幅均为光学遥感影像；第12到第18幅、第23幅以及第27到第34幅图像为SAR图像，其中第28幅是Radarsat图像，第23和31幅是SIR-C/X图像，第34幅是ERS-2图像；第22幅是IRS全色影像，第24幅为福卫二号卫星全色影像。
所用实验图像数据大小分布在200×200至1600×1600像素之间。

[0063] 将迭代阈值法、最大类间方差法、KSW熵又称一维Kapur熵方法、本发明方法分别与人工分割的水域参考图进行比较，给出了正确率、错误率、Kappa系数及运行时间等计算结果。其中人工分割的水域参考图是三位研究者分别进行分割的结果综合图。

[0064] Kappa系数是遥感分类图与参考图之间一致性或精度的度量，这种度量是通过主对角线和行列总数给出的概率一致性来表达的，它计算的是一个统计量K，是Kappa的一个估计值。将所述几种方法的分类图与人工参考图相比，如果K＞0.8说明比较方法和参考图的结果一致性很大或者精度很高；如果0.4＜K＜0.8表示一致性中等；如果K＜0.4表示一致性很差。

[0065] 正确检测数是参考图与比较方法中都被划分为水域的像元数；错误检测数是参考图划分为非水域而比较方法划分为水域的像元数；漏检数是参考图划分为水域而比较方法划分为非水域的像元数。错误率为错检率与漏检率之和；错检率为错误检测数除以参考图水域像元总数；漏检率为漏检数除以漏检数与比较方法水域像元总数之和。

[0066] 2.实验结果

[0067] 实验证明，本发明对SAR图像及满足微波水域反射特性的光学遥感影像均有效。附图3(a)和附图4(a)分别为表1中的第29和第11幅图，附图3(c)和附图4(c)分别为对应图像的水域识别的直观结果。表1中的35幅实验图像采用迭代阈值法、最大类间方差法及KSW熵方法同本发明方法的错误率、Kappa系数和本发明方法的运行时间等客观实验结果都列在表1中。从表1的相应列数据中可以看出，迭代阈值法、最大类间方差法及KSW熵方法的错误率远高于本发明方法，高出几倍甚至几十倍。35幅实验数据的Kappa的平均值数据充分说明了几种方法与参考图的一致程度。KSW熵方法平均Kappa系数最低只有
0.31，表示与参考图一致性很差，最大类间方差法的平均Kappa值为0.54，稍好于KSW熵方法，表示与参考图一致性稍差，迭代阈值法的平均Kappa值为0.69表示与参考图一致性较好，而本发明方法的平均Kappa值达到了0.91，这说明了本发明方法与参考图一致性最好，或者说本方法的精确度最高。从大量实验数据中可以看出，本发明方法比现有的方法识别效果更好，运行时间更快，检测精度更高这是其他现有方法难以达到的。

[0068] 表1实验数据结果比较

[0069]