CN114998724B - 一种抗云层干扰的湖泊蓝藻水华遥感监测方法 - Google Patents

Abstract

本发明公开了一种抗云层干扰的湖泊蓝藻水华遥感监测方法，包括以下步骤：1)构建水体像素、蓝藻水华像素和云像素样本集；2)计算水体像素、蓝藻水华像素和云像素的归一化系数；3)构建指标函数CBI；4)计算蓝藻水华像素阈值；5)提取蓝藻水华像素。本发明能够解决研究区域遥感影像中云层像素误判为蓝藻水华的问题，并能显著降低云层的误判率，从而提高遥感影像识别蓝藻水华的抗云层干扰能力。

Description

一种抗云层干扰的湖泊蓝藻水华遥感监测方法

技术领域

本发明涉及湖泊蓝藻水华遥感监测领域，具体涉及一种抗云层干扰的湖泊蓝藻水华遥感监测方法。

背景技术

河湖蓝藻在我国正呈蔓延趋势，成为主要污染来源。蓝藻快速生长聚集形成水华，导致水体生物多样性急剧下降，破坏水体景观和生态系统平衡，给周边经济带来了巨大的障碍，严重影响到区域环境，并制约当地的经济和社会的可持续发展，开展水华预防、治理研究应当从蓝藻监测开始。因此，快速、准确的掌握蓝藻水华分布信息，对水华预防、治理研究显得尤为重要。

目前，应用遥感影像提取蓝藻水华主要涉及到提取指标的问题，常见的提取指标主要有归一化植被指数(NDVI)、浮游藻类指数(FAI)、虚拟植被指数(VB-FAH)，虽然提取方法多样，但是原理基本一致，都是通过不同波段之间的组合，突出蓝藻水华的信息，然后通过设置一定的阈值将蓝藻水华提取出来。以上三种指标在一定程度上均可以达到蓝藻水华提取的目的，但是它们存在一定的问题，首先是每一种提取指标需要应用大气校正后的地表反射率数据，将对应的波段反射率数据带入到计算公式中即可突出遥感影像中蓝藻水华的信息。然而，不同的大气校正方法处理后的遥感影像，指标的计算结果存在显著差异，从而导致提取的蓝藻面积大相径庭。其次，云在光学遥感卫星中十分常见，而现有的蓝藻水华识别方法往往将一部分的云像素误判为蓝藻水华，降低了水华监测精度。

选择适合于特定研究对象的大气校正方法，是现有算法的必要前提，而不恰当的大气校正方法会进入较为严重的误差，往往需要实地测量水体光谱反射率，以选取适合的大气校正方法，耗费时间和人力，不仅如此，现有方法涉及图像预处理(包括大气校正)——指标计算——阈值划分——蓝藻水华确定，整个流程的耗时主要集中在大气校正环节。

因此，如何解决传统提取指标由于大气校正的不同和云像素存在的误判情况而导致的提取结果的差异性问题，成为了我国内陆湖泊中蓝藻水华遥感识别的迫切任务。

发明内容

本发明的目的是针对传统蓝藻水华提取指标将云像素误判为蓝藻的问题，提出了一种抗云层干扰的湖泊蓝藻水华遥感监测方法，以期能显著降低对云像素的误判率，提高遥感影像监测蓝藻水华的抗云层干扰能力，从而提高多光谱遥感卫星在内陆湖泊图像提取中的应用价值。

本发明为解决技术问题采用如下技术方案：

本发明一种抗云层干扰的湖泊蓝藻水华遥感监测方法的特点是按如下步骤进行：

步骤1、构建湖泊蓝藻水华提取指标；

步骤1.1构建水体、蓝藻水华、云像素样本集；

从研究区域的Landsat8 OL遥感影像中提取水体像素样本集合蓝藻水华像素样本集合/>和云像素样本集合/>其中，/>分别为第n个水体像素样本在红、绿、蓝和近红外波段的灰度值，/>分别为第n′个蓝藻水华像素样本在红、绿、蓝和近红外波段的灰度值，分别为第n″个云像素样本在红、绿、蓝和近红外波段的灰度值，N_wat、N_cyb和N_cld分别为集合中水体像素、蓝藻水华像素和云像素的个数，1≤n≤N_wat，1≤n′≤N_cyb，1≤n″≤N_cld；

利用卫星传感器A获取研究区域的另一遥感影像，分别提取所述水体像素样本集合蓝藻水华像素样本集合和云像素样本集合/>的同名点，得到水体像素同名点样本集合/>蓝藻水华像素同名点样本集合/>和云像素同名点样本集合其中，/>分别为第n个水体像素同名点样本在红、绿、蓝和近红外波段的灰度值，分别为第n′个蓝藻水华像素同名点样本在红、绿、蓝和近红外波段的灰度值，/>分别为第n″个云像素同名点样本在红、绿、蓝和近红外波段的灰度值；

步骤1.2计算湿度分量；

构建波段湿度分量系数矩阵{α^R、α^G、α^B、α^NIR}，其中，α^R、α^G、α^B和α^NIR分别为红、绿、蓝和近红外波段的湿度分量系数；

根据所述水体像素样本集合与所述波段湿度分量系数矩阵{α^R、α^G、α^B、α^NIR}，利用式(1)得到水体像素样本的湿度分量集合其中，y(n)为所述水体像素样本集合第n个水体像素样本的湿度分量；

步骤1.3、计算水体像素的归一化系数：

以所述水体像素同名点样本集合为观测点，利用式(2)所示的四元一次线性模型对湿度分量集合Y进行线性回归，得到湿度分量拟合函数y；

式(2)中，为水体像素在红波段的系数，/>为水体像素样本集合在红波段的灰度值，/>为水体像素在绿波段的系数，/>为水体像素样本集合在绿波段的灰度值，/>为水体像素在蓝波段的系数，/>为水体像素样本集合在蓝波段的灰度值，/>为水体像素在近红外波段的系数，/>为水体像素样本集合在近红外波段的灰度值；

令水体像素系数集合为

计算归一化的水体像素系数集合其中，A_wat表示水体像素系数集合的2范数，且/> 表示水体像素在红波段的归一化系数，/>表示水体像素在绿波段的归一化系数，/>表示水体像素在蓝波段的归一化系数，/>表示水体像素在近红外波段的归一化系数；

步骤1.4、计算云像素的归一化系数：

利用式(3)得到水体像素和云像素的相关系数γ_cld.wat：

式(3)中，表示云像素样本在红波段的灰度平均值，且/> 表示水体像素样本在红波段的灰度平均值，且/> 表示云像素样本在绿波段的灰度平均值，且/> 表示水体像素样本在绿波段的灰度平均值，且 表示云像素样本在蓝波段的灰度平均值，且/> 表示水体像素样本在蓝波段的灰度平均值，且/> 表示云像素样本在近红外波段的灰度平均值，且/> 表示水体像素样本在近红外波段的灰度平均值，且/>

利用式(4)-式(7)得到云像素系数集合

式(4)-式(7)中，为云像素在红波段的系数，/>为云像素在绿波段的系数，/>为云像素在蓝波段的系数，/>为云像素在近红外波段的系数；

计算归一化的云像素系数集合其中，A_cld表示云像素系数集合的2范数，且/> 表示云像素在红波段的归一化系数，/>表示云像素在绿波段的归一化系数，/>表示云像素在蓝波段的归一化系数，/>表示云像素在近红外波段的归一化系数；

步骤1.5、计算蓝藻水华像素的归一化系数：

利用式(8)得到蓝藻水华像素样本集合与水体像素样本集合的相关系数γ_cyb.wat：

式(8)中，表示蓝藻水华像素在红波段的灰度平均值，/> 表示蓝藻水华像素在绿波段的灰度平均值，/> 表示蓝藻水华像素在蓝波段的灰度平均值，/> 表示蓝藻水华像素在近红外波段的灰度平均值，

利用式(9)得到蓝藻水华像素与云像素的相关系数γ_cyb.cld：

利用式(10)-式(13)得到蓝藻水华像素系数集合

式(10)-式(13)中，为蓝藻水华像素在红波段的系数，/>为蓝藻水华像素在绿波段的系数，/>为蓝藻水华像素在蓝波段的系数，/>为蓝藻水华像素在近红外波段的系数；

计算归一化的蓝藻水华像素系数集合其中，A_cyb表示蓝藻水华像素系数集合的2范数，/> 表示蓝藻水华像素在红波段的系数，/>表示蓝藻水华像素在绿波段的系数，/>表示蓝藻水华像素在蓝波段的系数，/>表示蓝藻水华像素在近红外波段的系数；

步骤1.6、计算指标系数：

利用式(14)-式(17)得到在红、绿、蓝和近红外波段的指标系数k^R、k^G、k^B和k^NIR：

步骤1.7、利用式(18)构建指标函数CBI：

CBI＝k^R·λ^R+k^G·λ^G+k^B·λ^B+k^NIR·λ^NIR (18)

式(18)中，λ^R、λ^G、λ^B和λ^NIR分别为像素在红、绿、蓝和近红外波段的灰度值；

步骤2、确定阈值：

步骤2.1、利用式(18)分别计算所述水体像素同名点样本集合的指标函数CBI值集合/>蓝藻水华像素同名点样本集合/>的指标函数CBI值集合和云像素同名点样本集合/>的指标函数CBI值集合/>其中，CBI_wat(n)表示第n个水体像素同名点样本的指标函数CBI值，且/>CBI_cyb(n′)表示第n′个蓝藻水华像素同名点样本的指标函数CBI值，且CBI_cld(n″)表示第n″个云像素同名点样本的指标函数CBI值，且

步骤2.2、利用式(19)计算蓝藻水华像素阈值T：

式(18)中，max<·〉、min<·>分别取最大值和取最小值；

步骤3、蓝藻水华像素提取：

步骤3.1、构建研究区域的多波段遥感影像其中，p_mn表示第m行n列像素，且/> 和/>分别为第m行n列像素在红、绿、蓝和近红外波段的灰度值；

步骤3.2、利用式(18)计算所述研究区域的多波段遥感影像P中各像素的指标函数CBI值，构成指标函数CBI值集合其中，c_mn表示第m行n列像素p_mn的指标函数CBI值，且/>

步骤3.3、提取所述指标函数CBI值集合CBI^P中大于蓝藻水华像素阈值T的所有元素，构成指标函数CBI值集合其中，第i个元素/>即为指标函数CBI值集合CBI^P中的第m_i行n_i列的元素，从而得到所述研究区域多波段遥感影像P中的蓝藻水华像素集合/>其中，第i个元素/>即为所述研究区域多波段遥感影像P中的第m_i行n_i列像素，且1≤m_i≤m_I和1≤n_i≤n_I，1≤m_I≤M，1≤n_I≤N，m_I表示指标函数CBI值集合/>的总行数，n_I表示指标函数CBI值集合/>的总列数。

本发明与现有技术相比，本发明的有益效果在于：

1、云在光学卫星遥感影像中十分常见，而现有方法往往将一些云像素误判为蓝藻水华，严重降低水华监测精度。而本发明提出的蓝藻水华识别方法，通过Landsat8 OLI遥感影像数据的湿度分量系数，构建蓝藻水华提取指标，能够显著降低对云像素的误判率，从而提高了水华监测精度。

2、本发明提出的蓝藻水华识别方法，无需对遥感影像进行大气校正处理，使得整个流程的耗时显著降低，是一种快速的蓝藻水华识别方案。

附图说明

图1为本发明一种抗云层干扰的湖泊蓝藻水华遥感监测方法的流程图；

图2为本发明一种蓝藻水华监测方法所提取的蓝藻水华遥感监测结果图。

具体实施方式

本实施例中，一种抗云层干扰的湖泊蓝藻水华遥感监测方法，如图1所示，是按如下步骤进行：

步骤1、构建湖泊蓝藻水华提取指标；

步骤1.1构建水体、蓝藻水华、云像素样本集；

从研究区域的Landsat8 OL遥感影像中提取水体像素样本集合蓝藻水华像素样本集合/>和云像素样本集合/>其中，/>分别为第n个水体像素样本在红、绿、蓝和近红外波段的灰度值，/>分别为第n′个蓝藻水华像素样本在红、绿、蓝和近红外波段的灰度值，分别为第n″个云像素样本在红、绿、蓝和近红外波段的灰度值，N_wat、N_cyb和N_cld分别为集合中水体像素、蓝藻水华像素和云像素的个数，1≤n≤N_wat，1≤n′≤N_cyb，1≤n″≤N_cld；

利用卫星传感器A获取研究区域的另一遥感影像，分别提取所述水体像素样本集合蓝藻水华像素样本集合和云像素样本集合的同名点，得到水体像素同名点样本集合蓝藻水华像素同名点样本集合和云像素同名点样本集合其中，/>分别为第n个水体像素同名点样本在红、绿、蓝和近红外波段的灰度值，分别为第n′个蓝藻水华像素同名点样本在红、绿、蓝和近红外波段的灰度值，/>分别为第n″个云像素同名点样本在红、绿、蓝和近红外波段的灰度值；

步骤1.2计算湿度分量；

参考Chao Chen,Jiaoqi Fu,Shuai Zhang,Xin Zhao,Coastline informationextraction based on the tasseled cap transformation of Landsat-8 OLIimages.Estuarine,Coastal and Shelf Science,2019,217,281-291中的方法来构建波段湿度分量系数矩阵{α^R、α^G、α^B、α^NIR}，其中，α^R、α^G、α^B和α^NIR分别为红、绿、蓝和近红外波段的湿度分量系数；

根据体像素样本集合与所述波段湿度分量系数矩阵{α^R、α^G、α^B、α^NIR}，利用式(1)得到水体像素样本的湿度分量集合/>其中，y(n)为所述水体像素样本集合第n个水体像素样本的湿度分量；

步骤1.3、计算水体像素的归一化系数：

以所述水体像素同名点样本集合为观测点，利用式(2)所示的四元一次线性模型对湿度分量集合Y进行线性回归，得到湿度分量拟合函数y；

式(2)中，为水体像素在红波段的系数，/>为水体像素样本集合在红波段的灰度值，/>为水体像素在绿波段的系数，/>为水体像素样本集合在绿波段的灰度值，/>为水体像素在蓝波段的系数，/>为水体像素样本集合在蓝波段的灰度值，/>为水体像素在近红外波段的系数，/>为水体像素样本集合在近红外波段的灰度值；

令水体像素系数集合为

计算归一化的水体像素系数集合其中，A_wat表示水体像素系数集合的2范数，且/> 表示水体像素在红波段的归一化系数，/>表示水体像素在绿波段的归一化系数，/>表示水体像素在蓝波段的归一化系数，/>表示水体像素在近红外波段的归一化系数；

步骤1.4、计算云像素的归一化系数：

利用式(3)得到水体像素和云像素的相关系数γ_cld.wat：

式(3)中，表示云像素样本在红波段的灰度平均值，且/> 表示水体像素样本在红波段的灰度平均值，且/> 表示云像素样本在绿波段的灰度平均值，且/> 表示水体像素样本在绿波段的灰度平均值，且 表示云像素样本在蓝波段的灰度平均值，且/> 表示水体像素样本在蓝波段的灰度平均值，且/> 表示云像素样本在近红外波段的灰度平均值，且/> 表示水体像素样本在近红外波段的灰度平均值，且/>

利用式(4)-式(7)得到云像素系数集合

式(4)-式(7)中，为云像素在红波段的系数，/>为云像素在绿波段的系数，/>为云像素在蓝波段的系数，/>为云像素在近红外波段的系数；

计算归一化的云像素系数集合其中，A_cld表示云像素系数集合的2范数，且/> 表示云像素在红波段的归一化系数，/>表示云像素在绿波段的归一化系数，/>表示云像素在蓝波段的归一化系数，/>表示云像素在近红外波段的归一化系数；

步骤1.5、计算蓝藻水华像素的归一化系数：

利用式(8)得到蓝藻水华像素样本集合与水体像素样本集合的相关系数γ_cyb.wat：

式(8)中，表示蓝藻水华像素在红波段的灰度平均值，/> 表示蓝藻水华像素在绿波段的灰度平均值，/> 表示蓝藻水华像素在蓝波段的灰度平均值，/> 表示蓝藻水华像素在近红外波段的灰度平均值，

利用式(9)得到蓝藻水华像素与云像素的相关系数γ_cyb.cld：

利用式(10)-式(13)得到蓝藻水华像素系数集合

式(10)-式(13)中，为蓝藻水华像素在红波段的系数，/>为蓝藻水华像素在绿波段的系数，/>为蓝藻水华像素在蓝波段的系数，/>为蓝藻水华像素在近红外波段的系数；

计算归一化的蓝藻水华像素系数集合其中，A_cyb表示蓝藻水华像素系数集合的2范数，/> 表示蓝藻水华像素在红波段的系数，/>表示蓝藻水华像素在绿波段的系数，/>表示蓝藻水华像素在蓝波段的系数，/>表示蓝藻水华像素在近红外波段的系数；

步骤1.6、计算指标系数：

利用式(14)-式(17)得到在红、绿、蓝和近红外波段的指标系数k^R、k^G、k^B和k^NIR：

步骤1.7、利用式(18)构建指标函数CBI：

CBI＝k^R·λ^R+k^G·λ^G+k^B·λ^B+k^NIR·λ^NIR (18)

式(18)中，λ^R、λ^G、λ^B和λ^NIR分别为像素在红、绿、蓝和近红外波段的灰度值；

步骤2、确定阈值：

步骤2.1、利用式(18)分别计算所述水体像素同名点样本集合的指标函数CBI值集合/>蓝藻水华像素同名点样本集合/>的指标函数CBI值集合和云像素同名点样本集合/>的指标函数CBI值集合/>其中，CBI_wat(n)表示第n个水体像素同名点样本的指标函数CBI值，且/>CBI_cyb(n′)表示第n′个蓝藻水华像素同名点样本的指标函数CBI值，且CBI_cld(n″)表示第n″个云像素同名点样本的指标函数CBI值，且

步骤2.2、利用式(19)计算蓝藻水华像素阈值T：

式(18)中，max<·>、min<·>分别取最大值和取最小值，本实例中，计算的到的蓝藻水华提取阈值T为252；

步骤3、蓝藻水华像素提取：

步骤3.1、构建研究区域的多波段遥感影像其中，p_mn表示第m行n列像素，且/> 和/>分别为第m行n列像素在红、绿、蓝和近红外波段的灰度值；

步骤3.2、利用式(18)计算所述研究区域的多波段遥感影像P中各像素的指标函数CBI值，构成指标函数CBI值集合其中，c_mn表示第m行n列像素p_mn的指标函数CBI值，且/>

步骤3.3、提取所述指标函数CBI值集合CBI^P中大于蓝藻水华像素阈值T的所有元素，构成指标函数CBI值集合其中，第i个元素/>即为指标函数CBI值集合CBI^P中的第m_i行n_i列的元素，从而得到所述研究区域多波段遥感影像P中的蓝藻水华像素集合/>其中，第i个元素/>即为所述研究区域多波段遥感影像P中的第m_i行n_i列像素，且1≤m_i≤m_I和1≤n_i≤n_I，1≤m_I≤M，1≤n_I≤N，m_I表示指标函数CBI值集合/>的总行数，n_I表示指标函数CBI值集合/>的总列数；如图2所示，CBI值大于252的像素被分类为蓝藻水华，白色部分为云。/>

Claims (1)

1.一种抗云层干扰的湖泊蓝藻水华遥感监测方法，其特征在于，是按如下步骤进行：

步骤1、构建湖泊蓝藻水华提取指标；

步骤1.1构建水体、蓝藻水华、云像素样本集；

从研究区域的Landsat8 OL遥感影像中提取水体像素样本集合蓝藻水华像素样本集合/>和云像素样本集合/>其中，/>分别为第n个水体像素样本在红、绿、蓝和近红外波段的灰度值，/>分别为第n′个蓝藻水华像素样本在红、绿、蓝和近红外波段的灰度值，分别为第n″个云像素样本在红、绿、蓝和近红外波段的灰度值，N_wat、N_cyb和N_cld分别为集合中水体像素、蓝藻水华像素和云像素的个数，1≤n≤N_wat，1≤n′≤N_cyb，1≤n″≤N_cld；

利用卫星传感器A获取研究区域的另一遥感影像，分别提取所述水体像素样本集合蓝藻水华像素样本集合和云像素样本集合/>的同名点，得到水体像素同名点样本集合/>蓝藻水华像素同名点样本集合/>和云像素同名点样本集合其中，/>分别为第n个水体像素同名点样本在红、绿、蓝和近红外波段的灰度值，/>分别为第n′个蓝藻水华像素同名点样本在红、绿、蓝和近红外波段的灰度值，分别为第n″个云像素同名点样本在红、绿、蓝和近红外波段的灰度值；

步骤1.2计算湿度分量；

构建波段湿度分量系数矩阵{α^R、α^G、α^B、α^NIR}，其中，α^R、α^G、α^B和α^NIR分别为红、绿、蓝和近红外波段的湿度分量系数；

根据所述水体像素样本集合与所述波段湿度分量系数矩阵{α^R、α^G、α^B、α^NIR}，利用式(1)得到水体像素样本的湿度分量集合/>其中，y(n)为所述水体像素样本集合第n个水体像素样本的湿度分量；

步骤1.3、计算水体像素的归一化系数：

以所述水体像素同名点样本集合为观测点，利用式(2)所示的四元一次线性模型对湿度分量集合Y进行线性回归，得到湿度分量拟合函数y；

式(2)中，为水体像素在红波段的系数，/>为水体像素样本集合在红波段的灰度值，/>为水体像素在绿波段的系数，/>为水体像素样本集合在绿波段的灰度值，/>为水体像素在蓝波段的系数，/>为水体像素样本集合在蓝波段的灰度值，/>为水体像素在近红外波段的系数，/>为水体像素样本集合在近红外波段的灰度值；

令水体像素系数集合为

计算归一化的水体像素系数集合其中，A_wat表示水体像素系数集合的2范数，且/> 表示水体像素在红波段的归一化系数，/>表示水体像素在绿波段的归一化系数，/>表示水体像素在蓝波段的归一化系数，/>表示水体像素在近红外波段的归一化系数；

步骤1.4、计算云像素的归一化系数：

利用式(3)得到水体像素和云像素的相关系数γ_cld.wat：

式(3)中，表示云像素样本在红波段的灰度平均值，且/> 表示水体像素样本在红波段的灰度平均值，且/> 表示云像素样本在绿波段的灰度平均值，且/> 表示水体像素样本在绿波段的灰度平均值，且 表示云像素样本在蓝波段的灰度平均值，且/> 表示水体像素样本在蓝波段的灰度平均值，且/> 表示云像素样本在近红外波段的灰度平均值，且/> 表示水体像素样本在近红外波段的灰度平均值，且/>

利用式(4)-式(7)得到云像素系数集合

式(4)-式(7)中，为云像素在红波段的系数，/>为云像素在绿波段的系数，/>为云像素在蓝波段的系数，/>为云像素在近红外波段的系数；

计算归一化的云像素系数集合其中，A_cld表示云像素系数集合的2范数，且/> 表示云像素在红波段的归一化系数，/>表示云像素在绿波段的归一化系数，/>表示云像素在蓝波段的归一化系数，/>表示云像素在近红外波段的归一化系数；

步骤1.5、计算蓝藻水华像素的归一化系数：

利用式(8)得到蓝藻水华像素样本集合与水体像素样本集合的相关系数γ_cyb.wat：

式(8)中，表示蓝藻水华像素在红波段的灰度平均值，/> 表示蓝藻水华像素在绿波段的灰度平均值，/> 表示蓝藻水华像素在蓝波段的灰度平均值，/> 表示蓝藻水华像素在近红外波段的灰度平均值，

利用式(9)得到蓝藻水华像素与云像素的相关系数γ_cyb.cld：

利用式(10)-式(13)得到蓝藻水华像素系数集合

式(10)-式(13)中，为蓝藻水华像素在红波段的系数，/>为蓝藻水华像素在绿波段的系数，/>为蓝藻水华像素在蓝波段的系数，/>为蓝藻水华像素在近红外波段的系数；

计算归一化的蓝藻水华像素系数集合其中，A_cyb表示蓝藻水华像素系数集合的2范数，/> 表示蓝藻水华像素在红波段的系数，/>表示蓝藻水华像素在绿波段的系数，/>表示蓝藻水华像素在蓝波段的系数，/>表示蓝藻水华像素在近红外波段的系数；

步骤1.6、计算指标系数：

利用式(14)-式(17)得到在红、绿、蓝和近红外波段的指标系数k^R、k^G、k^B和k^NIR：

步骤1.7、利用式(18)构建指标函数CBI：

CBI＝k^R·λ^R+k^G·λ^G+k^B·λ^B+k^NIR·λ^NIR (18)

式(18)中，λ^R、λ^G、λ^B和λ^NIR分别为像素在红、绿、蓝和近红外波段的灰度值；

步骤2、确定阈值：

步骤2.1、利用式(18)分别计算所述水体像素同名点样本集合的指标函数CBI值集合/>蓝藻水华像素同名点样本集合/>的指标函数CBI值集合和云像素同名点样本集合/>的指标函数CBI值集合/>其中，CBI_wat(n)表示第n个水体像素同名点样本的指标函数CBI值，且/>CBI_cyb(n′)表示第n′个蓝藻水华像素同名点样本的指标函数CBI值，且CBI_cld(n″)表示第n″个云像素同名点样本的指标函数CBI值，且

步骤2.2、利用式(19)计算蓝藻水华像素阈值T：

式(18)中，max<·>、min<·>分别取最大值和取最小值；

步骤3、蓝藻水华像素提取：

步骤3.1、构建研究区域的多波段遥感影像其中，p_mn表示第m行n列像素，且/> 和/>分别为第m行n列像素在红、绿、蓝和近红外波段的灰度值；

步骤3.2、利用式(18)计算所述研究区域的多波段遥感影像P中各像素的指标函数CBI值，构成指标函数CBI值集合其中，c_mn表示第m行n列像素p_mn的指标函数CBI值，且/>

步骤3.3、提取所述指标函数CBI值集合CBI^P中大于蓝藻水华像素阈值T的所有元素，构成指标函数CBI值集合其中，第i个元素/>即为指标函数CBI值集合CBI^P中的第m_i行n_i列的元素，从而得到所述研究区域多波段遥感影像P中的蓝藻水华像素集合/>其中，第i个元素/>即为所述研究区域多波段遥感影像P中的第m_i行n_i列像素，且1≤m_i≤m_I和1≤n_i≤n_I，1≤m_I≤M，1≤n_I≤N，m_I表示指标函数CBI值集合/>的总行数，n_I表示指标函数CBI值集合/>的总列数。