CN109374537B - 城市黑臭水体识别方法及装置 - Google Patents

Abstract

本申请实施例公开了一种城市黑臭水体识别方法及装置，利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算水面待识别点在绿光波段的插值反射率；根据水面待识别点的绿光波段的遥感反射率和计算得到的插值反射率，以及红光波段的遥感反射率，计算黑臭水体分级指数；根据该黑臭水体分级指数判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级。不仅有效识别出黑臭水体，还实现了对黑臭水体的分级识别。

Description

城市黑臭水体识别方法及装置

技术领域

本申请涉及遥感影像处理技术领域，更具体地说，涉及一种城市黑臭水体识别方法及装置。

背景技术

城市建成区内受到污染达到一定程度的水体称为黑臭水体。城市黑臭水体不仅给群众带来了极差的感官体验，也是直接影响群众的生产生活。因此，需要对城市黑臭水体进行监测以便进行治理。

随着GF系列、ZY系列等国产高分辨率卫星的发射，我国遥感监测技术已达到亚米级水平，可以利用遥感手段大范围、实时、快速的监测城市黑臭水体。

然而目前的黑臭水体识别方法只能区分是否是黑臭水体，而对于黑臭水体，却不能区分是重度黑臭水体还是轻度黑臭水体。

发明内容

本申请的目的是提供一种城市黑臭水体识别方法及装置，以至少部分的克服现有技术中存在的技术问题。

为实现上述目的，本申请提供了如下技术方案：

一种城市黑臭水体识别方法，包括：

获取水面待识别点的遥感反射率；

利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率；

根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数；

根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级。

上述方法，优选的，所述根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数，包括：

将所述绿光波段的遥感反射率与所述插值反射率的差值，除以所述红光波段的遥感反射率，得到所述黑臭水体分级指数。

上述方法，优选的，所述根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级，包括：

若所述黑臭水体分级指数小于第一预设阈值，确定所述水面待识别点为重度黑臭水体；

若所述黑臭水体分级指数大于或等于所述第一预设阈值，且小于或等于第二预设阈值，确定所述水面待识别点为轻度黑臭水体；所述第二预设阈值大于所述第一预设阈值；

若所述黑臭水体分级指数大于所述第二预设阈值，确定所述水面待识别点为非黑臭水体。

上述方法，优选的，所述利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率，包括：采用线性内插方式计算所述插值反射率：

R'_rs,G＝R_rs,B+(R_rs,R-R_rs,B)*(λ_G-λ_B)/(λ_R-λ_B)

其中，R'_rs,G为所述绿光波段的插值反射率；R_rs,B为所述蓝光波段的遥感反射率；R_rs,R为所述红光波段的遥感反射率；λ_R表示红光波段的中心波长；λ_G表示绿光波段的中心波长；λ_B表示蓝光波段的中心波长。

上述方法，优选的，所述获取水面待识别点的遥感反射率，包括：

获取被监测地理区域的全色影像和多光谱影像；

对所述全色影像和所述多光谱影像进行融合，得到融合影像；

对所述融合影像进行正射校正，得到第一校正影像；

对所述第一校正影像进行第一大气校正，得到第二校正影像；

在所述第二校正影像中，分别利用水体指数法和近红外波段阈值法提取水体区域，将两种方法提取的水体区域求交集，得到目标水体区域；

在所述目标水体区域中获取水面待识别点的遥感反射率。

一种城市黑臭水体识别装置，包括：

获取模块，用于获取水面待识别点的遥感反射率；

第一计算模块，用于利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率；

第二计算模块，用于根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数；

分级模块，用于根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级。

上述装置，优选的，所述第二计算模块具体用于：将所述绿光波段的遥感反射率与所述插值反射率的差值，除以所述红光波段的遥感反射率，得到所述黑臭水体分级指数。

上述装置，优选的，所述分级模块具体用于：

若所述黑臭水体分级指数小于第一预设阈值，确定所述水面待识别点为重度黑臭水体；

若所述黑臭水体分级指数大于或等于所述第一预设阈值，且小于或等于第二预设阈值，确定所述水面待识别点为轻度黑臭水体；所述第二预设阈值大于所述第一预设阈值；

若所述黑臭水体分级指数大于所述第二预设阈值，确定所述水面待识别点为非黑臭水体。

上述装置，优选的，所述第一计算模块具体用于：采用线性内插方式计算所述插值反射率：

R'_rs,G＝R_rs,B+(R_rs,R-R_rs,B)*(λ_G-λ_B)/(λ_R-λ_B)

其中，R'_rs,G为所述绿光波段的插值反射率；R_rs,B为所述蓝光波段的遥感反射率；R_rs,R为所述红光波段的遥感反射率；λ_R表示红光波段的中心波长；λ_G表示绿光波段的中心波长；λ_B表示蓝光波段的中心波长。

上述装置，优选的，所述获取模块包括：

获取单元，用于获取被监测地理区域的全色影像和多光谱影像；

融合单元，用于对所述全色影像和所述多光谱影像进行融合，得到融合影像；

正射校正单元，用于对所述融合影像进行正射校正，得到第一校正影像；

大气校正单元，用于对所述第一校正影像进行第一大气校正，得到第二校正影像；

提取单元，用于在所述第二校正影像中，分别利用水体指数法和近红外波段阈值法提取水体区域，将两种方法提取的水体区域求交集，得到目标水体区域；

读取单元，用于在所述目标水体区域中读取水面待识别点的遥感反射率。

通过以上方案可知，本申请提供的一种城市黑臭水体识别方法及装置，利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算水面待识别点在绿光波段的插值反射率；根据水面待识别点的绿光波段的遥感反射率和计算得到的插值反射率，以及红光波段的遥感反射率，计算黑臭水体分级指数；根据该黑臭水体分级指数判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级。不仅有效识别出黑臭水体，还实现了对黑臭水体的分级识别。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例提供的城市黑臭水体识别方法的一种实现流程图；

图2为本申请实施例提供的获取水面待识别点的遥感反射率的一种实现流程图；

图3为本申请实施例提供的城市黑臭水体识别装置的一种结构示意图；

图4a是基于本申请提供的BOCI模型确定的样本与BOCI值的对应关系示意图；

图4b是基于现有技术提供的BOI模型确定的样本与BOI值的对应关系示意图；

图4c是基于现有技术提供的红绿波段比值法模型确定的样本与红绿波段比值法模型值的对应关系示意图；

图5a是基于本申请提供的BOCI模型确定的样本与BOCI值的对应关系的另一示意图；

图5b是基于现有技术提供的BOI模型确定的样本与BOI值的对应关系的另一示意图；

图5c是基于现有技术提供的红绿波段比值法模型确定的样本与红绿波段比值法模型值的对应关系的另一示意图。

说明书和权利要求书及上述附图中的术语“第一”、“第二”、“第三”“第四”等(如果存在)是用于区别类似的部分，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本申请的实施例能够以除了在这里图示的以外的顺序实施。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

发明人在实现本申请的过程中发现，黑臭水体遥感反射率在可见光蓝、绿、红这三个波段的变化相比一般水体更平缓，且黑臭水体在绿波段的基线高(绿波段的基线高是指：绿光波段的遥感反射率与绿光波段的插值反射率的差值，绿光波段的插值反射率的计算在后文中详细介绍)相比一般水体明显要小。利用绿波段的基线高可以将黑臭水体和一般水体区分开来，但这并不能区分重度黑臭水体与轻度黑臭水体。

发明人在实现本申请的过程中还发现，重度黑臭水体Ⅱ在红光波段的遥感反射率最高，轻度黑臭水体在红光波段的遥感反射率次之，一般水体在红光波段的遥感反射率低于轻度黑臭水体在红光波段的遥感反射率，重度黑臭水体Ⅰ在红光波段的反射率最低，但重度黑臭水体在可见光蓝、绿、红这三个波段的遥感反射率几乎相同(在光谱曲线上表现为几乎呈一条直线)，且绿光波段的基线高几乎为零。重度黑臭水体Ⅰ和重度黑臭水体Ⅱ统称为重度黑臭水体。发明人根据可见光波段的反射率将重度黑臭水体分为两类：将可见光波段反射率极低(即低于第一反射率阈值)的重度黑臭水体定义为重度黑臭水体Ⅰ，而将可见光波段反射率极高(即高于第二反射率阈值，第二反射率阈值高于第一反射率阈值)的重度黑臭水体定位为重度黑臭水体Ⅱ。

基于上述发现，本申请提取一种不仅可以快速识别黑臭水体，还可以定量评估黑臭水体污染程度的黑臭水体分级指数(Black and Odorous water ClassificationIndex，BOCI)模型。

请参阅图1，图1为本申请实施例提供的城市黑臭水体识别方法的一种实现流程图，可以包括：

步骤S11：获取水面待识别点的遥感反射率。

水面待识别点的遥感反射率可以直接通过遥感影像获取；或者，

水面待识别点的遥感反射率为：由水面待识别点的实测遥感反射率转换得到的卫星波段等效反射率。每一个水面采样点都可以作为一个水面待识别点。

步骤S12：利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算水面待识别点在绿光波段的插值反射率。

可以利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率，以及红光波段、绿光波段和蓝光波段的中心波长进行线性插值计算，得到水面待识别点在绿光波段的插值反射率。

步骤S13：根据水面待识别点的绿光波段的遥感反射率和插值反射率，以及水面待识别点的红光波段的遥感反射率，计算黑臭水体分级指数。

本申请实施例中，利用水面待识别点的绿光波段的实际遥感反射率，以及前述计算得到的绿光波段的插值反射率，计算水面待识别点的用于对水体进行分级识别的BOCI。

步骤S14：根据黑臭水体分级指数判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级。

本申请实施例中，黑臭水体分为两个等级：重度黑臭水体和轻度黑臭水体。可以结合不同的阈值判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级。

本申请提供的城市黑臭水体识别方法，利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算水面待识别点在绿光波段的插值反射率；根据水面待识别点的绿光波段的遥感反射率和计算得到的插值反射率，以及红光波段的遥感反射率，计算水面待识别点的黑臭水体分级指数；根据该黑臭水体分级指数判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级。不仅有效识别出黑臭水体，还实现了对黑臭水体的分级识别。

在一可选的实施例中，上述根据水面待识别点的绿光波段的遥感反射率和插值反射率，以及水面待识别点的红光波段的遥感反射率，计算黑臭水体分级指数的一种实现方式可以为：

将水面待识别点的绿光波段的遥感反射率与上述计算得到的插值反射率的差值，除以水面待识别点的红光波段的遥感反射率，得到黑臭水体分级指数BOCI。具体用公式可以表示为:

其中，R_rs,G为绿光波段的遥感反射率；R'_rs,G为绿光波段的插值反射率；R_rs,R为红光波段的遥感反射率。

通过该方法计算得到的BOCI，不仅缩小了重度黑臭水体Ⅰ与重度黑臭水体Ⅱ之间的光谱差异，还增大了重度黑臭水体、轻度黑臭水体、一般水体之间的光谱差异。

在一可选的实施例中，上述根据黑臭水体分级指数判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级的一种实现方式可以为：

若黑臭水体分级指数小于第一预设阈值，确定水面待识别点为重度黑臭水体。

若黑臭水体分级指数大于或等于上述第一预设阈值，且小于或等于第二预设阈值，确定水面待识别点为轻度黑臭水体；其中，第二预设阈值大于第一预设阈值。

若黑臭水体分级指数大于上述第二预设阈值，确定水面待识别点为非黑臭水体，即一般水体。

其中，第一预设阈值可以为0.12，第二预设阈值可以为0.26。

城市黑臭水体一般为悬浮物主导的黑臭水体，除了悬浮物主导的黑臭水体外，城市黑臭水体可能还包括其它类型的黑臭水体，如叶绿素主导的黑臭水体。上述阈值应用于悬浮物主导的黑臭水体的识别时，识别率最高，应用于其它类型的黑臭水体的识别时，识别率稍微低一些，但影响不大。基于此，可以适当调整阈值的大小，以便适应不同类型的黑臭水体。

在一可选的实施例中，利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算水面待识别点在绿光波段的插值反射率的一种实现方式可以为：

采用线性内插方式计算插值反射率：

R'_rs,G＝R_rs,B+(R_rs,R-R_rs,B)*(λ_G-λ_B)/(λ_R-λ_B)

其中，R'_rs,G为计算得到的绿光波段的插值反射率；R_rs,B为蓝光波段的遥感反射率；R_rs,R为红光波段的遥感反射率；λ_R表示红光波段的中心波长；λ_G表示绿光波段的中心波长；λ_B表示蓝光波段的中心波长。

若获取水面待识别点的遥感反射率所使用的遥感影像为GF2号卫星影像，则λ_R＝660nm，λ_G＝555nm，λ_B＝485nm。

在一可选的实施例中，上述获取水面待识别点的遥感反射率的一种实现流程图如图2所示，可以包括：

步骤S21：获取被监测地理区域的全色影像(Panchromatic image，PAN)和多光谱影像(Multi-spectral image，MS)。

待识别水域位于上述被监测地理区域内。全色影像和多光谱影像是同一时刻采集的被监测地理区域的高空间分辨率遥感影像。

步骤S22：对获取的全色影像和多光谱影像进行融合，得到融合影像。

可以利用ENVI软件自带的Gram-schmidt Pan Sharpening方法对多光谱影像和全色影像进行融合，该方法可以较好的保证融合前后影像的光谱形状和数值的一致性。具体在计算时，由于全色影像和多光谱影像均为DN(Digital Number)值影像(即影像中像元的值为DN值)，因此，使用全色影像和多光谱影像的DN值进行融合计算。计算得到融合影像也为DN值影像。

步骤S23：对融合影像进行正射校正，得到第一校正影像。

可以利用数字高程模型(Digital Elevation Model，DEM)辅助的区域网平差技术，对融合影像的成像几何模型进行基于有理函数模型的像方仿射变换修正，实现融合影像的正射校正并保证影像间几何位置的一致性。或者，

可以利用ENVI软件自带的正射校正模块对融合影像进行正射校正。

步骤S24：对第一校正影像进行第一大气校正，得到第二校正影像。

可以利用ENVI软件自带的大气校正模块对第一校正影像进行第一大气校正。

或者，

为了提高水体区域的提取精度，可以利用哨兵2A卫星拍摄的被监测地理区域的影像对第一校正影像进行第一大气校正。具体可以为：

对哨兵2A卫星拍摄的被监测地理区域的影像进行第二大气校正，得到第三校正影像。将该第三校正影像作为对第一校正影像进行大气校正的参考影像。可以利用欧空局提供的插件Sen2Cor对哨兵2A卫星拍摄的被监测地理区域的影像进行第二大气校正。

利用第三校正影像对第一校正影像进行第一大气校正，具体可以为：

步骤a、对第一校正影像进行重采样，得到分辨率与第三校正影像的分辨率相同的影像。

哨兵2A卫星拍摄的影像的空间分辨率通常较低，而第一校正影像是高空间分辨率影像。本申请实施例中，对第一校正影像进行重采样，以得到低空间分辨率影像，且该低空间分辨率影像的空间分辨率与第三校正影像的分辨率相同。

步骤b：分别对重采样得到的影像，以及第三校正影像进行不变地物搜索。

可以利用多元变化检测(Multivariate Alteration Detection，MAD)算法搜索重采样得到的影像的不变地物，以及利用多元变化检测算法搜索第三校正影像的不变地物。

步骤c：根据重采样得到的影像中，以及第三校正影像中的相同不变地物的数据，拟合线性方程。

重采样得到的影像为DN值影像，而大气校正后的哨兵2A卫星影像为反射率影像(即影像中像元的值为反射率)。因此，在拟合线性方程时，可以利用重采样得到的影像中不变地物的各个像元的DN值，以及第三校正影像中的相同不变地物的各个像元的反射率值，拟合线性方程。该拟合线性方程即是对第一校正影像进行第一大气校正的大气校正模型。

步骤S24：将第一校正影像的数据带入拟合线性方程，得到第二校正影像。

将第一校正影像中各个像元的DN值带入拟合线性方程，得到第二校正影像，该第二校正影像为反射率影像。

步骤S25：在第二校正影像中，分别利用NDWI水体指数法和近红外波段阈值法提取水体区域，将两种方法提取的水体区域求交集，得到目标水体区域。

可以利用ENVI软件里的band math工具，分别利用上述两种方法提取水体区域。具体可以为：

在band math工具中输入NDWI水体指数公式，并指定波段及对应的阈值，由bandmath工具根据上述NDWI水体指数公式，指定的波段及阈值对第二校正影像进行处理，得到与水体指数法对应的水体区域；

在band math工具中输入近红外波段计算公式，并指定近红外波段及对应的阈值，由band math工具根据上述近红外波段计算公式，指定的近红外波段及阈值对上述第二校正影像进行处理，得到与近红外波段阈值法对应的水体区域。

进一步的，为了提高水体区域的精度，还可以手动将求交集得到的水体区域中明显误判的像元(如，裸土，桥梁等)删除，将较为细小的河流进行补全，得到最终完整的河流水系区域，即目标水体区域。

步骤S26：在目标水体区域中读取水面待识别点的遥感反射率。

具体的，目标水体区域中，每一个像元都可以作为一个水面待识别点。

与方法实施例相对应，本申请还提供一种城市黑臭水体识别装置。本申请提供的城市黑臭水体识别装置的一种结构示意图如图3所示，可以包括：

获取模块31，第一计算模块32，第二计算模块33和分级模块34；其中，获取模块31用于获取水面待识别点的遥感反射率；

第一计算模块32用于利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率；

第二计算模块33用于根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数；

分级模块34用于根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级。

本申请提供的城市黑臭水体识别装置，利用水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算水面待识别点在绿光波段的插值反射率；根据水面待识别点的绿光波段的遥感反射率和计算得到的插值反射率，以及红光波段的遥感反射率，计算黑臭水体分级指数；根据该黑臭水体分级指数判断水面待识别点是否是黑臭水体，以及当水面待识别点是黑臭水体时，黑臭水体的等级。不仅有效识别出黑臭水体，还实现了对黑臭水体的分级识别。

在一可选的实施例中，第二计算模块33具体可以用于：将所述绿光波段的遥感反射率与所述插值反射率的差值，除以所述红光波段的遥感反射率，得到所述黑臭水体分级指数。

在一可选的实施例中，分级模块34具体可以用于：

若所述黑臭水体分级指数小于第一预设阈值，确定所述水面待识别点为重度黑臭水体；

若所述黑臭水体分级指数大于或等于所述第一预设阈值，且小于或等于第二预设阈值，确定所述水面待识别点为轻度黑臭水体；所述第二预设阈值大于所述第一预设阈值；

若所述黑臭水体分级指数大于所述第二预设阈值，确定所述水面待识别点为非黑臭水体。

在一可选的实施例中，第一计算模块32具体可以用于：采用线性内插方式计算所述插值反射率：

R'_rs,G＝R_rs,B+(R_rs,R-R_rs,B)*(λ_G-λ_B)/(λ_R-λ_B)

其中，R'_rs,G为所述绿光波段的插值反射率；R_rs,B为所述蓝光波段的遥感反射率；R_rs,R为所述红光波段的遥感反射率；λ_R表示红光波段的中心波长；λ_G表示绿光波段的中心波长；λ_B表示蓝光波段的中心波长。

在一可选的实施例中，获取模块31可以包括：

获取单元，用于获取被监测地理区域的全色影像和多光谱影像；

融合单元，用于对所述全色影像和所述多光谱影像进行融合，得到融合影像；

正射校正单元，用于对所述融合影像进行正射校正，得到第一校正影像；

大气校正单元，用于对所述第一校正影像进行第一大气校正，得到第二校正影像；

提取单元，用于在所述第二校正影像中，分别利用水体指数法和近红外波段阈值法提取水体区域，将两种方法提取的水体区域求交集，得到目标水体区域；

读取单元，用于在所述目标水体区域中读取水面待识别点的遥感反射率。

下面以2015年7月16-18日、2016年9月19-20日、2016年10月9-11日和2017年10月24-25日在沈阳市建成区范围内开展的四次野外水面试验所采集的光谱数据为例，将本申请的技术方案与现有技术的方案进行对比说明。本次对比使用的数据包括：在满堂河、辉山明渠、细河以及北运河和微山湖路附近判别和采集的69个黑臭水体样点的数据，在浑河和蒲河采集的47个一般水体样点的数据，共69+47＝116个采样点。

首先，从116个采样点中随机选择2/3样本(重度黑臭：32个；轻度黑臭：15个；一般水体：30个)进行阈值选取。将选择的样本数据带入黑臭水体识别模型得到阈值，具体的：

将选择的样本数据带入本申请的BOCI模型，如图4a所示，为基于BOCI模型确定的样本与BOCI值的对应关系示意图。本申请的BOCI模型识别一般水体与黑臭水体的阈值设为0.26，识别重度黑臭水体与轻度黑臭水体的阈值为0.12，本申请提供的BOCI模型不仅可以识别黑臭水体，还可以将重度黑臭水体与轻度黑臭水体区分开来。

将选择的样本数据带入姚月等提出的BOI模型(详见《基于高分二号的沈阳市黑臭水体遥感识别》)，如图4b所示，为基于BOI模型确定的样本与BOI值的对应关系示意图。文章给定的BOI阈值为0.065，认为小于0.065的为黑臭水体，大于0.065的为一般水体，但从图4b中可以看出，除了有与BOCI相同的误判点位之外，还有个别点位误判，根据样本数据重新修改阈值为0.075。

将选择的样本数据带入温爽等提出的红绿波段比值法模型(详见《基于高分影像的城市黑臭水体遥感识别:以南京为例》)，如图4c所示，为基于红绿波段比值法模型确定的样本与红绿波段比值法模型值的对应关系示意图。文章中给定的阈值为0.06，认为小于0.06的为黑臭水体，大于0.06的为一般水体，但从图4c中可以看出，0.06这个阈值无法有效的区分沈阳城市黑臭水体与一般水体，根据样本数据重新修改阈值为0.115。

通过对比分析可知，BOI模型和红绿波段比值法模型可以用于区分黑臭水体与一般水体，但是不能将轻度黑臭水体与重度黑臭水体区分开。而使用BOCI模型，不仅可以有效的识别黑臭水体，还可以将轻度黑臭水体与重度黑臭水体区分开。

下面将116个采样点中剩余的1/3样本(重度黑臭：15个；轻度黑臭：7个；一般水体：17个)用于精度验证。

将上述1/3样本数据带入BOCI模型，如图5a所述，为基于BOCI模型确定的样本与BOCI值的对应关系的另一示意图。阈值采用0.26，黑臭水体识别正确率为97.43％。

将选择的样本数据带入BOI模型，如图5b所示，为基于BOI模型确定的样本与BOI值的对应关系的另一示意图。阈值采用0.075，黑臭水体识别正确率为92.13％。

将选择的样本数据带入红绿波段比值法模型，如图5c所示，为基于红绿波段比值法模型确定的样本与红绿波段比值法模型值的对应关系的另示意图。阈值采用0.115，黑臭水体识别正确率为89.74％。

发明人研究发现，对于蓝光波段遥感反射率较高，绿-红波段遥感反射率斜率较大的情况，BOI算法与红绿波段比值法均不适用，而BOCI模型能够较好的识别。通过对比可以看出：BOCI模型识别黑臭水体的精度高于BOI模型和红绿波段比值法模型。

由于BOCI模型对黑臭水体进行了轻重分级，下面基于混淆矩阵对BOCI模型进行分级精度评价，如表1(表1中数据基于图5a得到)所示，总体精度为87.18％，Kappa为0.80。

表1基于混淆矩阵的BOCI模型精度评价

表1中，每一行代表了预测类别(即BOCI模型的识别结果)，每一行的总数表示预测为该类别的数据的数目；每一列代表了数据的真实归属类别(即调研结果)，每一列的数据总数表示该类别的数据实际的数目。每一行中的数值表示真实数据被预测为该类的数目。例如，第一行第一列的17表示有17个实际为一般水体的样本被预测为一般水体。第一行第三列的1表示有1个实际为重度黑臭水体的样本被错误预测为一般水体。第二行第三列的3表示有3个实际为重度黑臭水体的样本被错误的预测为轻度黑臭水体。

生产者精度：是指BOCI模型将所有样本中正确分为A类的样本数(对角线值)与A类真实参考总数(混淆矩阵中A类列的总和)的比率。如本例中，轻度黑臭水体有7个真实样本，其中6个正确分类，因此轻度黑臭水体的生产者精度是6/7＝85.71％。

用户精度：是指正确分到A类的样本总数(对角线值)与BOCI模型将所有样本的样本分为A类的样本总数(混淆矩阵中A类行的总和)比率。如本例中，轻度黑臭水体有6个正确分类，总共划分为轻度黑臭水体的有9个，所以轻度黑臭水体的用户精度是6/9＝66.67％。

总体精度为：被正确分类的样本总和除以总样本数。被正确分类的样本数目沿着混淆矩阵的对角线分布，总样本数等于所有真实样本总数。如本例中，正确分类的样本总和为：17+6+11＝34，总样本数为17+7+15＝39，总体精度为：34/39＝87.18％。

Kappa为：它是通过把所有真实样本总数(39)乘以混淆矩阵对角线数据的和(34)，再减去某一类样本总数与被分成该类样本总数之积对所有类别求和的结果(17*18+7*9+15*12＝549)，再除以样本总数的平方(39*39)减去某一类样本总数与被分成该类样本总数之积对所有类别求和的结果(549)所得到的，即Kappa＝(39*34-549)/(39*39-549)＝0.80。

Kappa值越高，说明BOCI的结果与实际调研结果的一致性越高。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

应当理解，本申请实施例中，从权、各个实施例、特征可以互相组合结合，都能实现解决前述技术问题。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (6)

1.一种城市黑臭水体识别方法，其特征在于，包括：

获取水面待识别点的遥感反射率；

利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率；

根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数；

根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级；

所述利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率，包括：采用线性内插方式计算所述插值反射率：

R'_rs,G＝R_rs,B+(R_rs,R-R_rs,B)*(λ_G-λ_B)/(λ_R-λ_B)

其中，R'_rs,G为所述绿光波段的插值反射率；R_rs,B为所述蓝光波段的遥感反射率；R_rs,R为所述红光波段的遥感反射率；λ_R表示红光波段的中心波长；λ_G表示绿光波段的中心波长；λ_B表示蓝光波段的中心波长；

所述根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数，包括：

将所述绿光波段的遥感反射率与所述插值反射率的差值，除以所述红光波段的遥感反射率，得到所述黑臭水体分级指数。

2.根据权利要求1所述的方法，其特征在于，所述根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级，包括：

若所述黑臭水体分级指数小于第一预设阈值，确定所述水面待识别点为重度黑臭水体；

若所述黑臭水体分级指数大于或等于所述第一预设阈值，且小于或等于第二预设阈值，确定所述水面待识别点为轻度黑臭水体；所述第二预设阈值大于所述第一预设阈值；

若所述黑臭水体分级指数大于所述第二预设阈值，确定所述水面待识别点为非黑臭水体。

3.根据权利要求1所述的方法，其特征在于，所述获取水面待识别点的遥感反射率，包括：

获取被监测地理区域的全色影像和多光谱影像；

对所述全色影像和所述多光谱影像进行融合，得到融合影像；

对所述融合影像进行正射校正，得到第一校正影像；

对所述第一校正影像进行第一大气校正，得到第二校正影像；

在所述第二校正影像中，分别利用水体指数法和近红外波段阈值法提取水体区域，将两种方法提取的水体区域求交集，得到目标水体区域；

在所述目标水体区域中获取水面待识别点的遥感反射率。

4.一种城市黑臭水体识别装置，其特征在于，包括：

获取模块，用于获取水面待识别点的遥感反射率；

第一计算模块，用于利用所述水面待识别点的蓝光波段的遥感反射率和红光波段的遥感反射率计算所述水面待识别点在绿光波段的插值反射率；

第二计算模块，用于根据所述水面待识别点的绿光波段的遥感反射率和插值反射率，以及所述红光波段的遥感反射率，计算黑臭水体分级指数；

分级模块，用于根据所述黑臭水体分级指数判断所述水面待识别点是否是黑臭水体，以及当所述水面待识别点是黑臭水体时，黑臭水体的等级；

所述第一计算模块具体用于：采用线性内插方式计算所述插值反射率：

R'_rs,G＝R_rs,B+(R_rs,R-R_rs,B)*(λ_G-λ_B)/(λ_R-λ_B)

其中，R'_rs,G为所述绿光波段的插值反射率；R_rs,B为所述蓝光波段的遥感反射率；R_rs,R为所述红光波段的遥感反射率；λ_R表示红光波段的中心波长；λ_G表示绿光波段的中心波长；λ_B表示蓝光波段的中心波长

所述第二计算模块具体用于：将所述绿光波段的遥感反射率与所述插值反射率的差值，除以所述红光波段的遥感反射率，得到所述黑臭水体分级指数。

5.根据权利要求4所述的装置，其特征在于，所述分级模块具体用于：

若所述黑臭水体分级指数小于第一预设阈值，确定所述水面待识别点为重度黑臭水体；

若所述黑臭水体分级指数大于或等于所述第一预设阈值，且小于或等于第二预设阈值，确定所述水面待识别点为轻度黑臭水体；所述第二预设阈值大于所述第一预设阈值；

若所述黑臭水体分级指数大于所述第二预设阈值，确定所述水面待识别点为非黑臭水体。

6.根据权利要求4所述的装置，其特征在于，所述获取模块包括：

获取单元，用于获取被监测地理区域的全色影像和多光谱影像；

融合单元，用于对所述全色影像和所述多光谱影像进行融合，得到融合影像；

正射校正单元，用于对所述融合影像进行正射校正，得到第一校正影像；

大气校正单元，用于对所述第一校正影像进行第一大气校正，得到第二校正影像；

提取单元，用于在所述第二校正影像中，分别利用水体指数法和近红外波段阈值法提取水体区域，将两种方法提取的水体区域求交集，得到目标水体区域；

读取单元，用于在所述目标水体区域中读取水面待识别点的遥感反射率。

Images