基于序列非线性滤波的遥感影像水体专题信息提取方法

摘要

本发明属于遥感技术领域，具体涉及一种基于序列非线性滤波的遥感影像水体专题信息提取方法。具体步骤为：遥感影像的预处理与序列非线性滤波，滤波结果影像区域标定和各标定区域各项特征值的计算，水体目标区域选择准则的建立，水体区域的自动提取及基于形态学的水体专题提取信息后处理四个步骤，最后生成满足通用GIS软件要求的专题信息矢量数据。本发明从灰度数学形态学基本运算入手，构成一种序列非线性滤波模型，直接进行非线性滤波，提取遥感影像中的水体区域；根据水体的相关知识，对序列非线性滤波提取出的水体区域，计算不同区域的各项特征值根据水体成像特性和水体区域相关知识，建立水体目标区域选择准则，从而实现遥感影像上水体专题信息的自动快速提取，本发明受噪声和误差影响较少，鲁棒性强，处理提取水体过程全自动进行，无须人工干预，能处理的数据广泛，耗时少。

说明书

技术领域

本发明属于遥感技术领域，具体涉及一种基于序列非线性滤波的遥感影像水体专题信息提取方法。

背景技术

水体是重要的专题信息，它对于地学分析有重要的价值，是一类基础的专题信息数据。对于各地理信息系统的管理与建设有重要意义。利用遥感影像进行水体信息提取是一种快速、且现势性好的方式。目前利用遥感影像提取村镇专题信息的方法有：人工方法、半自动跟踪提取方法，以及自动模式分类识别方法；这些方法中人工方法、半自动获取方法需要大量的人工操作效率低；利用多光谱信息分类识别方法是一种重要的方法，但现有的分类识别方法不是专门针对遥感影像上的水体信息提取而研究的，而且这种识别分类方法对于一般分辨率(空间分辨率20m，或更低)的多波段遥感影像才比较适用；对于高分辨率(分辨率达到5m及以上)且波段较少(4个波段或更少)的遥感影像，其对水体的分类识别提取结果不能满足后续成图及应用要求。目前，在高分辨率遥感影像分类识别方面，德国的eCognition系统的分类方法是比较好的一种，其基于多尺度分割，分类时需要有分类协议，每一次分类的分类协议用于其它分类时需经过修改，最终实现分类过程的自动化，这里的分类协议类似于知识支持。其特点是分类的细化，但对分类协议有较大依赖(有协议库支持较好)，并且也不是专门针对水体这类信息而进行的，因此应用于水体信息提取时效率不高。为提高水体提取效率，避免提取过程中的人工干扰，实现一种全自动的水体提取方法，本方法针对上述的情况产生了利用灰度数学形态学运算，构成一种序列非线性滤波模型，直接进行非线性滤波，有效提取遥感影像中的水体区域；在此基础上，根据水体的相关知识，对上面提取出的水体区域，计算不同标定区域的各项特征值(面积，区域灰度均值，区域的灰度直方图)，根据水体成像特性和水体区域相关知识，建立水体目标区域选择准则，以此进行水体区域自动选择提取，得到水体专题信息。

发明内容

本发明的目的在于，针对现有的利用高分辨率遥感影像进行水体专题信息提取方法存在的受各种细节信息干扰严重、提取效率低、效果差等缺点，提出一种基于序列非线性滤波的遥感影像水体专题信息提取方法，本发明在保留利用遥感影像进行水体专题信息提取现势性强、重复性高等优点的同时，对高分辨率多波段遥感影像的水体信息提取更加适用，其效率更高。克服了一般提取方法需要先验知识，对噪音及细节干扰信息敏感的缺点，利用灰度数学形态学运算来构成一种序列非线性滤波模型，直接对高分辨率遥感影像进行非线性滤波，统计标定目标区域的特征值，根据水体所具有的特征，来进行遥感影像水体专题信息自动识别。

本发明提出的基于序列非线性滤波的遥感影像水体专题信息提取方法，通过从灰度数学形态学基本运算入手，构成一种序列非线性滤波模型，直接进行非线性滤波，提取遥感影像中的水体区域；在此基础上，根据水体的相关知识，对序列非线性滤波提取出的水体区域，计算不同区域的各项特征值(面积，区域灰度均值，区域的灰度直方图)，根据水体成像特性和水体区域相关知识，建立水体目标区域选择准则，从而实现遥感影像上水体专题信息的自动快速提取，具体包括：遥感影像的预处理与序列非线性滤波，滤波结果影像区域标定和各标定区域各项特征值的计算，水体目标区域选择准则的建立，水体区域的自动提取及基于形态学的水体专题提取信息后处理四个步骤，具体如下：

(1)遥感多波段影像的预处理与序列非线性滤波

根据水体光谱特性和遥感影像波段设置，选择红外或红波段、蓝或绿波段作为备选波段，根据公式(1)进行预处理，并对预处理结果做直方图均衡化，生成待处理数据图层，

BW₀＝B_L+k×B_L/(B_H+w)        (1)

其中，B_L表示红外或红波段亮度值，B_H表示蓝或绿波段亮度值，BW₀表示预处理后的像素亮度值，k为比例变换系数，w是为了避免除零问题，而对B_H加上的一个很小的正数。

采用序列非线性滤波方法，对生成待处理数据图层首先进行最大值滤波处理，消除影像中非连续面状分布的暗的像素，同时呈面状连续分布的水体区域也被削弱；然后进行中值滤波进行降噪处理；进而采用最小值滤波处理恢复最大值滤波处理中被削弱的水体区域。

(2)对滤波结果影像进行区域标定、计算各个标定区域的各项特征值

采用区域生长算法，对步骤(1)中的结果影像先标定候选水体像素组成连续区域，然后计算各个区域的各项统计特征：区域面积、区域灰度均值、区域灰度直方图，

区域灰度直方图为：

$\left(\begin{array}{cc}H\left[i\right]=\underset{R\left[i\right]}{\Sigma }1& i=\mathrm{0,1},...,255\end{array}\right)---\left(2\right)$

区域面积A为：

$A=\underset{i=0}{\overset{255}{\Sigma }}H\left[i\right]---\left(3\right)$

区域灰度均值：

$\mathrm{GM}=\frac{1}{A}\underset{i=1}{\overset{255}{\Sigma }}[H\left[i\right]*i]---\left(4\right)$

直方图峰值灰度级：

GP＝PI，H[PI]＝MAX{H[i]，i＝0，1，...，255}        (5)

其中，R[i]为区域中灰度值为i的像素集合，H[i]为直方图第i级灰度的像素个数，A为区域面积，GM为区域灰度均值，GP为直方图峰值灰度级；

(3)水体目标区域选择准则的建立

①区域面积A<A0的区域为非水体区域；②区域灰度均值GM>GM0的区域为非水体；③区域峰值灰度GP>区域均值灰度GM的区域为非水体区域。其中A0及GM0取值均有较大的冗余度，可取A0＝100，GM0＝64。

(4)水体区域的自动提取，基于形态学的水体专题提取信息后处理

根据步骤(3)建立的水体区域准则，及各个区域的特征值，对标记的各个区域进行判定，从而自动剔除杂波及其他目标，得到水体目标区域；通过形态开、闭运算，对生成的专题图像进行滤波，除去比结构元素小的特定图像细节，在形态学滤波的基础上对提取的矢量进行长度及面积统计，对提取出的矢量长度及所围面积小于阈值的进行剔除。阈值可调整，对于高分辨率航空遥感影像一般长度阈值为25，面积阈值为500，最后提取的专题信息分层表示为ArcGIS软件的形文件，即shp文件格式。

以下对本发明方法作进一步的限定，具体如下：

1、遥感多波段影像的预处理与序列非线性滤波

根据水体光谱特性和遥感影像波段设置，选择红外(IR)及蓝(B)两个波段作为备选波段；针对不同遥感影像，如果没有波段设置或缺少上述波段，也可选择红(R)代替红外(IR)或绿(G)代替蓝(B)波段。则可根据下式进行预处理并对预处理结果作直方图均衡化，生成待处理数据图层。

BW₀＝B_L+k×B_L/(B_H+w)

其中，B_L表示红外(IR)或红(R)波段亮度值，B_H表示蓝(B)或绿(B)波段亮度值，BW₀表示预处理后的像素亮度值，k为比例变换系数，w是为了避免除零问题，而对B_H加上的一个很小的正数。

应用灰度形态学算子与非线性统计滤波算子，在灰度形态学运算中定义一个特殊的结构元素，将灰度形态学运算转化为非线性滤波处理，通过序列非线性滤波，将形态学灰度膨胀运算转化为最大值滤波处理，灰度腐蚀运算转化为最小值滤波处理。采用组合处理步骤，对遥感影像首先进行最大值滤波处理，消除影像中非连续面状分布的暗的像素，同时呈面状连续分布的水体区域也被削弱；然后进行中值滤波进行降噪处理；进而采用最小值滤波处理恢复最大值滤波处理中被削弱的水体区域。处理过程中邻域大小可以灵活选取，也可以连续进行多次最大值滤波，最小值滤波时也连续进行相同次数的处理，则对于连续面状分布的水体目标可以得到保留，而消除掉大部分非水体区域，提取遥感影像中水体目标区域。

2、对滤波结果影像进行区域标定、计算各个标定区域的各项特征值

采用区域生长算法，对(1)中的结果影像先标定候选水体像素组成连续区域，然后计算各个区域的各项统计特征：区域面积、区域灰度均值、区域灰度直方图。以此为自动建立各个区域的判断规则做准备，避免了影像全局阈值设定困难的问题。

3、水体目标区域选择准则的建立

根据水体成像特性和水体区域相关知识，建立水体目标区域选择准则，根据水体为连续面状分布，建立区域面积判定规则，对面积小于额定像素的区域判定为非水体区域；根据水体区域亮度比较统一，其直方图应为单峰形状，建立灰度分布判定规则，判定出非水体区域。这里亮度判断可以设定比较宽松的标准，结合直方图单峰特征，先排除掉明显非水体区域。再根据区域直方图峰值灰度和区域均值灰度的关系，自动判断水体区域。自动确定以下水体选择准则：1)区域面积A<A0的区域为非水体区域；2)区域灰度均值GM>GM0的区域为非水体；3)区域峰值灰度GP>区域均值灰度GM的区域为非水体区域。其中A0及GM0取值均有较大的冗余度，可取A0＝100，GM0＝64。

4、水体区域的自动提取，基于形态学的水体专题提取信息后处理

根据以上建立的水体区域准则，及各个区域的特征值，对标记的各个区域进行判定，从而自动剔除杂波及其它目标，得到水体目标区域；在此基础上，为了消除结果上因细小地物以及噪声影响，而产生的细小地物图斑，通过形态开、闭运算，对生成的专题图像进行滤波(结构元素的大小可以选择有3×3、5×5，7×7，15×15)，除去比结构元素小的特定图像细节，同时保证不产生全局的集合失真，在形态学滤波的基础上对提取的矢量进行长度及面积统计，对提取出的矢量长度及所围面积小于阈值的进行剔除。阈值可调整(对于高分辨率航空遥感影像一般长度阈值为25，面积阈值为500)，最后提取的专题信息分层表示为ArcGIS软件的形文件，即shp文件格式。

与现有技术相比，本发明采用了序列非线性滤波方法，避免了传统方法二值化中的阈值计算确定，速度更快，而且对高分辨率卫星影像更加适用，效率更高，采用区域直方图形状及灰度峰值均值相对关系判断方法，有效的避免了水体判断时的阈值问题，避免提取过程中作业者的人工干预，具有更好的鲁棒性和自适应特性。本发明可以应用于高分辨率多波段遥感影像的水体专题信息快速自动提取，在水体专题信息获取与变化检测方面发挥作用。

附图说明

图1为本发明实施例1中TM(专题成像仪)图像数据除热红外波段以外的其它谱段1、2、3、4、5、7，其中：(a)为1波段谱段(蓝)，(b)为2波段谱段，(c)为3波段谱段，(d)为4波段谱段，(e)为5波段谱段，(f)为7波段谱段(红外)。

具体实施方式

以下结合一个实施例子对具体实现方法进行说明，即要对一幅高分辨率多波段遥感影像进行水体专题信息识别，获得满足ArcGIS软件要求的水体专题信息矢量数据(.shp文件)。

实施例1：

这里取某一TM(专题成像仪)数据——图象大小为2516×2395像素，含7个波段，其已具有地理编码信息——图象左上角点坐标为(285950，3475350)，像素分辨率为25米，图1显示的是除热红外波段以外的其它谱段1、2、3、4、5、7。由于数据量极大，例子中(1)、(2)、(3)中的显示数据是该图象局部自图象像素行列坐标(753，723)至(768，738)范围的16×16大小图像块的处理结果。

(1)遥感多波段影像的预处理和序列非线性滤波

根据水体光谱特性和遥感影像波段设置，选择红外(IR)及蓝(B)两个波段作为备选波段；针对不同遥感影像，如果没有波段设置缺少上述波段，也可选择红(R)代替红外(IR)或绿(G)代替蓝(B)波段。定义B_L表示上述红外(IR)或红(R)波段亮度值，B_H表示上述蓝(B)或绿(B)波段亮度值，BW₀表示预处理后的像素亮度值。则可以下式进行预处理，并对预处理结果作直方图均衡化生成待处理数据图层。

BW₀＝B_L+k*B_L/(B_H+w)

其中，参数k及w取值都有很大的冗余度，这里取k＝20，w＝0.1。如果BW₀>255，则限制BW₀＝255。

然后，采用序列非线性滤波方法，针对上述预处理后图层，可根据公式(6)～(8)依次进行具有非线性特性的最大值滤波，中值滤波及最小值滤波。

BW₁＝MAX{U₈(BW₀)}         (6)

BW₂＝MED{U₈(BW₁)}         (7)

BW₃＝MIN{U₈(BW₂)}         (8)

其中，U₈(·)表示像素的8邻域，MAX(·)，MED(·)，MIN(·)分别表示取邻域内像素的最大值、中值及最小值，BW₀，BW₁，BW₂，BW₃分别表示预处理后图层、最大值滤波后、中值滤波后及最小值滤波后的像素亮度值。

这里，选用TM图像的第1波段(蓝)和第7波段(红外)作为BH和BL波段灰度数据。表1及表2为选取的一个16×16大小区域影像波段数据示例。表3及表4为线性变换及直方图均衡化后结果数据。表5为根据公式(6)～(8)进行序列非线性滤波后结果数据，其中已将大于128的像素灰度置为背景灰度255。

表1蓝色波段

102  104  108  110  107  105  106  103  105  110  109  103  102  100 100  102

100  104  107  107  106  106  105  102  104  107  107  103  104  102 100  103

100  105  110  110  107  106  108  110  107  110  109  104  102  101 99   102

97   103  107  105  105  106  107  106  104  106  107  104  101  101 102  103

101  102  101  98   103  105  104  104  103  101  103  105  104  102 103  102

101  100  100  102  104  103  100  100  101  101  101  104  104  103 104  102

97   99   102  103  103  103  100  99   98   99   103  103  100  102 104  102

97   100  100  99   99   102  103  102  101  101  101  102  99   101 103  101

99   99   100  98   97   100  103  104  101  102  102  102  101  102 102  100

100  100  100  98   99   100  102  101  102  103  102  101  104  105 102  99

100  99   99   98   100  99   99   100  102  101  102  102  103  103 101  99

99   99   101  101  99   99   99   100  100  101  104  101  100  104 102  101

100  101  101  100  99   101  100  100  101  100  99   101  102  105 103  104

99   100  99   99   100  100  100  99   100  102  100  102  103  102 101  101

99   100  101  101  100  101  100  100  100  102  102  100  102  102 102  102

99   100  103  101  100  102  104  103  102  102  102  99   101  101 102  103

                    表2红外波段

22  33  53  60  61  53  46  51  59  57  50  33  25  27  30  26

33  39  59  68  64  64  54  42  49  49  60  43  20  33  29  21

32  42  61  68  57  48  53  52  48  51  60  38  15  25  15  12

29  38  51  60  54  49  52  54  51  56  57  33  13  15  9   8

31  34  27  30  41  49  41  34  30  30  36  37  25  17  14  9

24  26  24  13  21  30  27  25  23  19  25  39  33  16  14  22

10  13  18  14  12  9   12  17  20  24  21  21  21  11  8   17

5   3   7   16  19  15  10  18  25  25  12  3   3   4   3   4

3   3   0   4   15  20  19  26  26  14  5   0   0   2   1   0

2   4   2   3   7   16  21  26  16  4   2   5   2   2   1   1

2   3   3   1   0   5   13  13  6   1   2   4   0   2   3   0

1   3   4   0   1   3   3   3   2   3   1   0   1   1   1   2

2   1   3   2   1   4   2   0   2   2   2   2   3   3   1   2

0   0   2   1   0   0   1   1   2   1   2   2   1   2   2   0

0   3   3   0   1   2   1   3   4   0   1   1   1   2   2   1

3   3   0   0   3   5   11      3   2   2   0   3   4   1   3

表3线性转换后结果

26  39  62  70  72  63  54  60  70  67  59  39  29  32  35  31

39  46  70  80  76  76  64  50  58  58  71  51  23  39  34  25

38  49  72  80  67  57  62  61  56  60  70  45  17  29  18  14

34  45  60  71  64  58  61  64  60  66  67  39  15  17  10  9

37  40  32  36  48  58  48  40  35  35  42  44  29  20  16  10

28  31  28  15  25  35  32  29  27  22  29  46  39  19  16  26

12  15  21  16  14  10  14  20  24  28  25  25  25  13  9   20

6   3   8   19  22  17  11  21  29  29  14  3   3   4   3   4

3   3   0   4   18  23  22  30  31  16  5   0   0   2   1   0

2   4   2   3   8   19  25  31  19  4   2   5   2   2   1   1

2   3   3   1   0   6   15  15  7   1   2   4   0   2   3   0

1   3   4   0   1   3   3   3   2   3   1   0   1   1   1   2

2   1   3   2   1   4   2   0   2   2   2   2   3   3   1   2

0   0   2   1   0   0   1   1   2   1   2   2   1   2   2   0

0   3   3   0   1   2   1   3   4   0   1   1   1   2   2   1

3   3   0   0   3   5   1   1   3   2   2   0   3   4   1   3

表4线性转换并经直方图均衡化后结果

60  110  182  203  207  185  157  176  203  195  173  110  72   85   96   81

110 132  203  221  214  214  187  144  170  170  205  147  50   110  92   56

106 141  207  221  195  167  182  179  164  176  203  129  41   72   42   38

92  129  176  205  187  170  179  187  176  193  195  110  39   41   34   33

103 113  85   100  138  170  138  113  96   96   120  126  72   45   40   34

67  81   67   39   56   96   85   72   64   49   72   132  110  44   40   60

35  39   47   40   38   34   38   45   53   67   56   56   56   37   33   45

29  27   32   44   49   41   34   47   72   72   38   27   27   28   27   28

27  27   0    28   42   50   49   76   81   40   28   0    0    25   22   0

25  28   25   27   32   44   56   81   44   28   25   28   25   25   22   22

25  27   27   22   0    29   39   39   30   22   25   28   0    25   27   0

22  27   28   0    22   27   27   27   25   27   22   0    22   22   22   25

25  22   27   25   22   28   25   0    25   25   25   25   27   27   22   25

0   0    25   22   0    0    22   22   25   22   25   25   22   25   25   0

0   27   27   0    22   25   22   27   28   0    22   22   22   25   25   22

27  27   0    0    27   28   22   22   27   25   25   0    27   28   22   27

         表5序列非线性滤波后结果

255  255  255  255  255  255  255  255  255  255  255  255  255  255  255  255

255  255  255  255  255  255  255  255  255  255  255  255  255  255  255  255

255  255  255  255  255  255  255  255  255  255  255  255  255  255  255  255

255  255  255  255  255  255  255  255  255  255  255  255  255  255  255  33

255  255  255  255  255  255  255  255  255  255  255  255  255  255  255  34

255  255  255  255  255  255  255  255  255  255  255  255  255  255  40   60

255  255  255  255  255  255  255  255  255  255  255  255  255  255  33   45

29   27   32   44   49   255  255  255  255  255  255  255  27   28   27   28

27   27   0    28   42   50   255  255  255  255  28   0    0    25   22   0

25   28   25   27   32   44   255  255  255  28   25   28   25   25   22   22

25   27   27   22   0    29   255  255  30   22   25   28   0    25   27   0

22   27   28   0    22   27   27   27   25   27   22   0    22   22   22   25

25   22   27   25   22   28   25   0    25   25   25   25   27   27   22   25

0    0    25   22   0    0    22   22   25   22   25   25   22   25   25   0

0    27   27   0    22   25   22   27   28   0    22   22   22   25   25   22

27   27   0    0    27   28   22   22   27   25   25   0    27   28   22   27

(2)对滤波结果影像进行区域标定、计算各个标定区域的各项特征值

对非线性滤波结果影像采用区域生长算法进行区域标定，然后统计各个标定区域的灰度直方图，根据下面4式，分别计算各个区域的面积、灰度均值和直方图峰值灰度。

$\left(\begin{array}{cc}H\left[i\right]=\underset{R\left[i\right]}{\Sigma }1& i=\mathrm{0,1},...,255\end{array}\right)$

$A=\underset{i=0}{\overset{255}{\Sigma }}H\left[i\right]$

$\mathrm{GM}=\frac{1}{A}\underset{i=1}{\overset{255}{\Sigma }}[H\left[i\right]*i]$

GP＝PI，H[PI]＝MAX{H[i]，i＝0，1，...，255}

其中，R[i]为区域中灰度值为i的像素集合，H[i]为直方图第i级灰度的像素个数，A为区域面积，GM为区域灰度均值，GP为直方图峰值灰度级。

(3)水体目标区域选择准则的建立

根据第(2)步计算的特征值，可以自动确定以下水体选择准则：1)区域面积A<A0的区域为非水体区域；2)区域灰度均值GM>GM0的区域为非水体；3)区域峰值灰度GP>区域均值灰度GM的区域为非水体区域。其中A0及GM0取值均有较大的冗余度，这里取A0＝100，GM0＝64。

(4)水体区域的自动提取，基于形态学的水体专题提取信息后处理

根据(3)建立的水体区域准则以及各个区域的特征值，自动对标记的各个区域进行判定，从而自动剔除杂波及其它目标，从备选区域中选择出水体区域。在此基础上，通过形态开、闭运算，对选出的水体区域进行形态学滤波(结构元素的大小可以选择有3×3、5×5，7×7，15×15)，去掉微小矢量分支(除去比结构元素小的特定图像细节)，同时保证不产生全局的集合失真，在形态学滤波的基础上，提取区域边缘矢量，对提取的矢量进行长度及面积统计，矢量长度及所围面积小于阈值的进行剔除。阈值可调整(对于高分辨率航空遥感影像一般长度阈值为25，面积阈值为500)，最后提取的专题信息分层表示为ArcGIS软件的形文件，即shp文件格式。表6为生成的某一地物专题矢量数据示例。

表6生成的某一地物专题矢量数据

点号专题矢量结点坐标(X，Y)

1~4   300600，3452250   300550，3452200   300475，3452200   300475，3452075

5~8   300450，3452050   300450，3451950   300425，3451925   300400，3451925

9~12    300350，3451875   300350，3451825   300325，3451800   300325，3451750

13~16   300350，3451725   300350，3451600   300325，3451575   300325，3451525

17~20   300350，3451500   300350，3451450   300375，3451425   300375，3451150

21~24   300350，3451125   300350，3451075   300375，3451050   300475，3451050

25~28   300500，3451075   300550，3451075   300575，3451100   300600，3451100

29~32   300625，3451125   301000，3451125   301025，3451150   301175，3451150

33~36   301200，3451125   301250，3451125   301250，3451225   301300，3451275

37~40   301375，3451275   301450，3451350   301575，3451350   301700，3451475

41~44   301900，3451475   301925，3451450   302150，3451450   302175，3451475

45~48   302175，3451550   302150，3451550   302100，3451600   302100，3451650

49~52   302075，3451675   302075，3451725   302100，3451750   302150，3451750

53~56   302150，3451850   302125，3451875   302075，3451875   302050，3451900

57~60   302000，3451900   301900，3452000   301900，3452025   301825，3452100

61~64   301800，3452100   301775，3452125   301750，3452125   301725，3452150

65~68   301700，3452150   301675，3452175   301600，3452175   301575，3452150

69~72   301550，3452150   301525，3452125   301500，3452125   301450，3452075

73~76   301425，3452075   301400，3452050   301375，3452050   301350，3452025

77~80   301300，3452025   301275，3452050   301250，3452050   301225，3452075

81~84   301025，3452075   301000，3452050   300850，3452050   300775，3452125

85~88   300775，3452150   300750，3452175   300750，3452250   300600，3452250