一种基于Sentinel-2影像的山区水体提取方法

一种基于sentinel-2影像的山区水体提取方法
技术领域
1.本发明属于遥感影像智能识别技术领域，尤其涉及一种基于sentinel-2影像的山区水体提取方法。


背景技术：

2.随着地球气候的变化，洪涝、水体滑坡等灾害频繁发生，对人民生命财产造成了重大损失，因此，快速准确地获取地表水体信息对洪涝、水体滑坡等灾害的监测具有重要意义。利用光学遥感影像获取水体信息是一种快速有效的方法，近年来，基光学遥感影像提取水体方法已由最初的人工目视解译发展到现在的自动分类提取水体方法，目前，常用的自动分类提取水体方法克服了人工目视解译不能大区域提取水体的缺点，该类方法主要包括基于像素分类的阈值法和基于目标分类的分类方法。阈值法更多的应用于中、高分辨率影像上，主要是基于地物光谱特征建立分类模型和水体指数，其算法简单有效，在湖泊等大面积水体提取中广泛应用，但不适用细小水体提取和山体水体提取。分类方法是基于影像光谱、纹理和空间特征来建立提取水体的方法，主要用于高分辨率影像中，常用的方法主要是svm方法、决策树分类法和面向对象方法，其具有提取更精细水体的优点，但难以将阴影和水体有效区分，易对提取山区水体造成混淆，并且该方法的算法复杂，需要建立训练样本，提取速度会因样本选择而降低。
3.水体指数法能实现山区水体的提取，目前常用的水体指数主要有：
①
归一化差异水体指数(ndwi)，能够最大程度抑制植被信息，突出水体信息，但是难以将建筑和土壤与水体区分，且受冰雪、云以及山体阴影的影响较大。
②
改进的归一化差异水体指数(mndwi)，能最大程度抑制建筑和土壤信息，突出水体特征，在城市水体提取中具有较好的效果。
③
增强水体指数(ewi)，能够抑制建筑、土壤和植被信息，突出半干旱区域的水体，但是易与阴影、雪混淆且与河滩难以区分，也不易提取山区水体。
④
阴影水体指数(swi)，具有区分水体和阴影的优点，同时也可以降低积雪和山体阴影的影响。然而，上述的水体指数方法在山区水体及细小水体提取中存在漏提现象，且存在将冰雪与云当作水体提取的不足，提取精度不高，在预测山体滑坡等自然灾害时存在误差。


技术实现要素：

4.针对上述背景技术中指出的不足，本发明提供了一种基于sentinel-2影像的山区水体提取方法，旨在解决上述背景技术中现有提取山区水体难度大，存在漏提现象，提取精度不高等的问题。
5.为实现上述目的，本发明采用的技术方案是：
6.一种基于sentinel-2影像的山区水体提取方法，方法如下：
7.步骤一：在sentinel-2msi影像上选取多种大小性质不同的个水体样本，统计分析多种水体样本在sentinel-2msi影像中的大气顶层反射率值，并计算各水体样本分别在不同影像波段的大气顶层反射率值；作为优选，所述影像波段包括b
11
(短波红外波段1)、b
12
(短波红外波段2)、b2(蓝光波段)、b3(绿光波段)、b4(红光波段)、b5(植被红边波段1)、b6(植被红边波段2)、b7(植被红边波段3)、b8(近红外波段)、b9(水蒸气波段)。
8.步骤二：根据各水体样本在影像不同波段的大气顶层反射率值，以与直角坐标系的横轴为45
°
夹角的直线为界限，将大气顶层反射率值分为两种分布趋势的情况，根据两种分布趋势将sentinel-2msi影像上的水体分为两类水体；然后选取可以反映两类水体的影像波段以及影像波段对应的大气顶层反射率值，将其反射率值和波段在新的直角坐标系下进行表示，并在新的直角坐标系中构建不同的三角形。
9.其中，三角形的构建方法具体为：以水体样本的影像波段对应的最大大气顶层反射率值和最小大气顶层反射率值的平均值为纵坐标，以影像波段为横坐标。
10.该步骤中，划分的两类水体包含的影像波段分别为：
11.第ⅰ类水体：b
12
、b3、b4、b5、b6；
12.第ⅱ类水体：b
12
、b2、b3、b9。
13.步骤三：根据两类水体各自构建的三角形特点，选择初始水体指数计算方法，并利用三角形相似定理，分别计算两类水体的水体指数，再将两类水体的水体指数求和，从而得到针对sentinel-2影像的水体指数，即三角形水体指数，利用所述三角形水体指数计算sentinel-2影像得到增强的水体图，然后通过直方图获取最佳阈值提取水体，获得最终的水体信息分布图。
14.优选地，上述划分的第ⅰ类水体计算水体指数时选择的初始水体指数计算公式如下：
15.twi1＝(b
3-b
12
)/(b3+b
12
)。
16.第ⅱ类水体计算水体指数时选择的初始水体指数计算公式如下：
17.twi2＝(b
9-b8)/(b8+b9)。
18.上述山区水体提取方法也可用于湖泊、水库、河流或者城市地表水的提取。
19.相比于现有技术的缺点和不足，本发明具有以下有益效果：
20.本发明无需建立训练样本，提取速度不受样本选择而变化，能够将山体阴影和水体有效区分。本发明在山区水体提取过程中避免了漏提现象，提取精度高，能够有效减少预测山体滑坡等自然灾害的误差率。
附图说明
21.图1是在sentinel-2msi影像上选取的11种水体样本及sentinel-2msi影像图。
22.图2是sentinel-2msi影像中水体样本的大气顶层反射率值统计分布图及各水体样本分别在不同影像波段的大气顶层反射率值统计分布图。
23.图3是两类水体的大气顶层反射率值与影像波段在直角坐标系中构建的三角形示意图。
24.图4是实施例1brahmaputra河不同区域处的sentinel-2msi影像图及进行twi指数法提取水体的结果图。
25.图5是实施例2和实施例3研究区的详细位置图。
26.图6是实施例2中青海省局部研究区域水体的不同水体指数提取结果图。
27.图7是实施例3中西藏局部研究区域的不同水体指数的结果图。
具体实施方式
28.为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
29.以sentinel-2影像作为本发明的遥感影像进行分析说明，在sentinel-2msi影像上选取11种大小性质不同的水体样本，包括河流、湖泊、水库和细小河流，如图1所示，统计分析11个水体样本在sentinel-2msi影像中的大气顶层反射率值，统计结果如图2(k)所示，可知不同性质的水体在各波段的反射率是存在相近但也表现出不同的特性。图2(k)中水体反射率走势表现出两种趋势，一类水体的反射率在波段b
11
到波段b2呈现下降趋势，在波段b8到波段b9表现出大幅度的增加趋势。而另一种水体的反射率在波段b
11
到波段b2呈现逐渐增加的趋势，在波段b8到波段b9没有明显的增加趋势且变化平稳。计算各水体样本分别在b
11
(短波红外波段1)、b
12
(短波红外波段2)、b2(蓝光波段)、b3(绿光波段)、b4(红光波段)、b5(植被红边波段1)、b6(植被红边波段2)、b7(植被红边波段3)、b8(近红外波段)、b9(水蒸气波段)的像素大气顶层反射率值，结果如图2(a)-(j)所示，图2中(a)、(b)、(c)、(d)、(e)、(f)、(g)、(h)、(i)、(j)分别为波段2、波段3、波段4、波段5、波段6、波段7、波段8、波段9、波段11、波段12上不同水体样本的反射率值。由图2(a)-(j)可以看出，水体样本的反射率值也表现为两种集中分布趋势，都沿直角坐标系右下方0点处至左上方的方向以45
°
为界限依据反射率强度将水体样本分为两类。
30.因此，根据直角坐标系下各水体样本在不同影像波段的像素大气顶层反射率值的分布情况，以与直角坐标系的横轴为45
°
夹角的直线，沿直角坐标系右下方0点处至左上方的方向，将sentinel-2msi影像的水体分为两类水体，两类水体的反射率值及对应的影像波段划分如表1所示：
31.表1sentinel-2a影像的水体
[0032][0033]
其中第ⅰ类水体的特点是在sentinel-2msi影像上颜色偏浅绿色，如样本water1-7所示；第ⅱ类水体特点是在sentinel-2msi影像上颜色偏暗，如样本water8-11所示。第ⅰ类水体在图2(k)中的水体反射率值在波段b3和波段b5有最高值，在b
12
和b6处有最低值；第ⅱ类水体在图2(k)中的水体反射率值在波段b
12
和波段b9反射强，在波段b2和b3处吸收强。根据以上特点分别建立两类水体在直角坐标系中的相应波段及大气顶层反射率值的三角形，如图3所示。图3(a)是第ⅰ类水体构建的三角形，(b)是第ⅱ类水体构建的三角形。其中，图3(a)的b5和b3是根据第一类水体样本的最大值和最小值的平均值；图3(b)的b9是根据第二类水体样本的最大值和最小值的平均值。如此，就可以将其转化为数学问题，利用数学特性建立水体指数。
[0034]
根据两类水体的三角形特点，选择两类水体的初始水体指数计算方法，并利用三角形相似定理，分别计算两类水体的水体指数，再将两类水体的水体指数求和，从而得到针
对sentinel-2影像的水体指数，即三角形水体指数(twi)，基于所得到的三角形水体指数确定sentinel-2影像上的山区水体。详细计算过程如下：
[0035]
第ⅰ类水体光谱特性符合改进的归一化水体指数(mndwi)研究方法，该方法最大程度抑制了建筑和土壤信息，突出水体特征，在城市水体提取中具有好的效果。在sentinel-2影像中mndwi计算公式如下：
[0036]
mndwi＝(b
3-b
12
)/(b3+b
12
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0037]
因此，本发明基于改进的归一化水体指数(mndwi)对构成该指数的(g-sir)利用三角形(图3(a)所示)原理进行了替换来计算第ⅰ类水体的水体指数。故计算第ⅰ类水体的水体指数时选择的初始水体指数计算公式如下：
[0038]
twi1＝(b
3-b
12
)/(b3+b
12
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0039]
图3(a)中，根据三角形相似判断定理可知δaco≌δdeo，δbeo≌δhfo，因此根据相似三角形相似定理：相似三角形的对应边成比例，可得到在δaco和δdeo中：
[0040][0041]
由公式(3)可得：
[0042][0043]
在δbeo和δhfo中：
[0044][0045]
由公式(5)可得：
[0046][0047]
将公式(4)代入公式(6)可得：
[0048][0049]
由于图3(a)中有：
[0050]b6-b5＝(b
12-b5)-(b
12-b6)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0051]
将公式(4)和公式(7)带入到公式(8)中可得到：
[0052][0053]
将公式(9)合并得：
[0054]
[0055]
由公式(10)可得：
[0056][0057]
因此，由twi1＝(b
3-b
12
)/(b3+b
12
)结合公式(11)可以改写成以下公式(12)，即第ⅰ类水体的水体指数：
[0058][0059]
分析第ⅱ类水体光谱反射率分布情况，根据植被指数的构建原理，基于sentinel-2影像利用b8(近红外波段)和b9(水蒸气波段)建立第ⅱ类水体指数(twi2)。构建的第ⅱ类水体的初始水体指数计算公式如下：
[0060]
twi2＝(b
9-b8)/(b8+b9)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0061]
由于波段b9在sentinel-2影像中是水蒸气波段，其不具有稳定性和实用性，所以利用三角形相似原理将b9用其他波段替换。由图3(b)中的δbgh可知：(b
9-b8)＝(b
9-b2)-(b
8-b2)，据此将公式(13)改为：
[0062][0063]
图3(b)中有δohb≌δodc和δaeg≌δbhg，由三角形相似定理可知：
[0064][0065]
通过公式(15)推导可得：
[0066][0067][0068]
由公式(15)和公式(17)可知：
[0069][0070]
将公式(18)变换可得到：
[0071][0072]
令|ae|＝a，|bh|＝c，公式(19)可改写成：
[0073][0074]
因此，可以得到b9为：
[0075][0076]
将公式(19)和(21)代入(14)式，即得第ⅱ类水体的水体指数：
[0077][0078]
依据上述计算过程分别求得了两类水体的水体指数计算公式，因此，总的三角形水体指数为两类水体指数求和的结果：
[0079][0080]
|de|、|be|、|hf|、|ac|均表示的是大小，为了便于使用，对系数进行了四舍五入。其中，|de|≈0.021，|be|≈0.041，|hf|≈0.031，|ac|≈0.032，a≈0.055，c≈0.199，基于水体样本反射率是小于1，因此，对取倒数变为取倒数变为则基于sentinel-2影像的三角形水体指数twi定义为：
[0081][0082]
在sentinel-2影像中，b
12
是短波红外波段，b2是蓝光波段，b3是绿光波段，b5是植被红边波段1，b6是植被红边波段2和b8是近红外波段。twi水体指数指在消除非水像素，能够实现在高海拔山区水体的提取，抑制山体阴影，雪等非水像素，该指数也可以实现抑制城市背景区域的深色建筑表面，提取城市水体。
[0083]
实施例1
[0084]
选取2021年7月28日的位于印度区域的brahmaputra河的sentinel-2a影像数据为实验点，brahmaputra河的北部区域海拔高且有较多的山地分布，山间有较小的湖泊和水洼，brahmaputra河的区域的具体位置如图4(a)和(d)所示，在gee平台采用twi指数方法计算图4(a)和(d)得到增强的水体分别如图4(b)和(e)所示，本研究采用otsu算法自动提取twi指数的阈值，获得最终水体信息分布分别如图(c)和(f)所示。对比水体提取前后图可知，利用本发明twi水体指数方法提取的水体信息分布位置和面积与sentinel-2msi影像展示的实际水体信息分布情况几乎一致，因此，twi水体指数方法具有良好的分辨效果，能够
很好的提取了山体间的水体，有效降低了噪声的影响，水体—非水体边界更加平滑，获得更加完整的河流、湖泊分布图。
[0085]
实施例2
[0086]
为定性评价本发明方法的有效性，体现本发明twi指数在水体信息提取方面的优势，选择当前经典的六种水体指数作为对比进行水体提取分析，经典的六种水体指数为：湿度指数(lswi)、改进的归一化差异水体指数(mndwi)、归一化差异水体指数(ndwi)、自动提取水体指数(awei)、阴影水体指数(swi)和增强水体指数(ewi)。选取2021年10月13日的西藏区域的sentinel-2a影像数据，西藏区域的水体包括河流、湖泊和其它常规水体，以及含有沉积物、盐和冰的特殊水体，且有大量冰雪存在，这些条件是检验twi指数有效性的最好条件。研究区的详细的位置如图5所示，西藏的研究区域位于东经89.88～91.76，北纬31.06～32.37。
[0087]
在gee平台上使用lswi、mndwi、ndwi、awei、swi、ewi和本分发明twi方法提取西藏研究区域的水体。各水体指数提取结果如图6所示。
[0088]
图6(a)为研究区域的谷歌地图，图6(b)为图6(a)的sentinel-2a影像数据，可以看出，影像数据所覆盖的地物有雪、湖水、山体、阴影和冰，经过ewi、awei、swi、lswi、mnwi、mndwi、twi水体指数获取的水体增强图分别如图6(c)、(d)、(e)、(f)、(g)、(h)、(i)所示。图6(d)、(e)和(f)显示水体与背景分离效果不明显，图6(c)、(g)和(h)水体增强效果明显优于图6(d)、(e)和(f)的水体增强性能。从增强的水体性能中可以看出，lswi、mndwi、ndwi、awei、swi和ewi难以识别甚至错过了一些结冰的水体(图中正方形区域)，并且细小河流识别不完全甚至未识别为水体(图中椭圆区域所示)，而twi解决了这些不足，可以较好的提取冰块水体，也能有效地将水体与雪区分以及增强了细小河流。由此可知，本发明twi可以有效实现山区湖水、冰冻水及河流的提取，且能够将雪、云和山体阴影与水体区分。
[0089]
实施例3
[0090]
选取2021年10月13日的青海省局部区域的sentinel-2a影像数据，研究区的详细的位置如图5所示，青海省的研究区域为东经97.46～98.38，北纬34.54～35.16的鄂灵湖、星星海和龙日阿错的区域。
[0091]
在gee平台上使用lswi、mndwi、ndwi、awei、swi、ewi和本分发明twi方法提取青海省局部研究区域的水体。各水体指数提取结果如图7所示。图7(a)和(b)分别为sentinel-2msi影像和sentinel-2msi影像上a区域的影像，图7(d)、(e)、(h)、(i)、(l)、(m)和(p)分别为sentinel-2msi影像的ewi、awei、swi、lswi、ndwi、mndwi和twi的提取结果图，图7(c)、(f)、(g)、(j)、(k)、(n)和(o)分别为sentinel-2msi影像上a区域影像的ewi、awei、swi、lswi、ndwi、mndwi和twi的提取结果图。
[0092]
由图7可知：各类水体指数都能增强水体，但是存在明显的差异，效果不同，图中黑色圆圈处可以看到swi(图(h))和lswi(图(i))水体与背景区分效果不佳，其次是awei(图(e))，与awei、swi和lswi相比，ewi(图(d))、mndwi(图(m))和ndwi(图(l))水体增强明显，能够分离湖水和背景，但是湖水面积有一些未曾增强。twi(图(p))可以有效地将水体特征与背景分离，且几乎可以将湖水面积完全增强。图(b)是细小的河流分布，综合比较图7中(c)、(f)、(g)、(j)、(k)、(n)和(o)，可知，ewi(图(c))、swi(图(g))、lswi(图(j))、mndwi(图(n))、ndwi(图(k))和awei(图(f))均不能有效增强山区细小河流，而利用twi(图(o))可以突出山
区河流，增强山区河流和背景。
[0093]
误差说明：
[0094]
上述选取西藏和青海省作为验证试验区，因google earth影像具有亚米级分辨率的优势，所以以google earth影像作为参考影像验证提取精度，在google earth影像上随机产生一定数量的采样点，对sentinel-2msi影像基于7种指数得到的水体分布进行精度评价，建立混淆矩阵，计算错漏分率、总体精度和kappa系数等定量评价指标，如表2和表3所示：
[0095]
表2水体指数评价指标
[0096][0097]
表3水体指数评价指标
[0098][0099][0100]
由表2和表3可以看出，以西藏和青海作为研究区域，相比较其它几种水体指数，twi水体识别结果的总体精度和kappa系数整体更优，总体精度为98％，错误识别和遗漏识别率低，分别是2.5％和2.3％，对山区水体提取识别效果较好。
[0101]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。