湖泊作为地表主要的生态组成部分，在调节局部气候、调蓄水资源和改善生态环境等方面有着重要作用^[1]，同时也为人类的生产生活提供保障。随着人类的活动日益加剧，水体富营养化问题已经严重威胁到湖泊生态安全^[2]。近年来我国主要的大型湖泊都出现过大面积藻华现象，如巢湖自2008年起一直保持富营养化水平^[3]，太湖在2004年后藻华面积快速增加，在2007年甚至一度达到湖面积的一半左右^[4]，滇池、洞庭湖等湖泊近十年均频繁出现藻华的现象^[5-6]。藻华暴发会对周围居民的生产生活和生态环境造成严重影响。因此，如何快速并准确获取藻华暴发的范围对改善湖泊周边生态环境，显得尤为重要^[7]。

随着遥感技术的不断发展，不同类型的光学卫星相继升空(如中国的高分系列卫星^[8]，欧空局Sentinel-2卫星^[9]，美国Landsat8、Landsat9卫星^[10]等)使水体环境的遥感监测得到长足发展，大范围、长时序、高精度的遥感监测内陆水体^[11]已成为可能。当藻华暴发时，暴发藻华的水体与清洁水体的光谱特征出现差异^[12]，藻华水体对红、蓝波段的吸收增加且对光的后向散射增强，尤其是近红外波段，会出现植被特有的“陡坡效应^[12]”现象，这种现象能够区分藻华水体和清洁水体，因此，“陡坡效应”可作为区分藻华水体和清洁水体的重要依据。目前水体表层藻华的遥感识别方法主要有：基于“陡坡效应”的单波段算法、近红外与可见光波段差异的差值法^[13]，比值法^[14]，归一化植被指数^[15](normalized difference vegetation index, NDVI)，增强型植被指数^[16](enchanted vegetation index, EVI)，浮游藻类指数^[17](floating algae index, FAI)等。

在利用遥感技术提取藻华的过程中，采用不同种类遥感影像提取藻华时，藻华光谱敏感指数阈值具有一定差异，甚至同一卫星不同时相的影像提取藻华时，藻华光谱敏感指数阈值也不相同^[18]，可见藻华光谱敏感指数提取藻华阈值存在不确定性^[19]，因此，对藻华的光谱敏感指数的阈值确定成为关键环节之一。针对该关键环节许多研究人员提出了各种不同的阈值确定方法。房旭等^[20]通过多源数据验证了高分辨率卫星对小型水库藻华的联合监测能力，在研究中采用NDVI指数，通过像元坡度计算确定NDVI提取藻华的阈值。Oyama等^[21]使用FAI指数法提取日本Nishiura、Kitaura和Hachiro三个小型湖泊的藻华情况，采用的VCL(visual cyanobacteria index，VCI)指标为确定FAI阈值范围，但是阈值具有一定的主观成分。Liu等^[22]以Landsat8和Sentinel-2为数据源，使用FAI指数提取云南省星云湖藻华范围，研究中利用FAI直方图确定阈值，但这种方法对于较小的藻华斑块提取结果并不理想。目前NDVI已经广泛应用于藻华的识别中，且采用坡度计算获取NDVI阈值的方法也得到验证。相对于NDVI，FAI对环境条件的改变敏感度较低，并且由于增加了短红外波段有效地降低了部分大气和薄云的影响，对藻华的识别有较高的稳定性。但如何确定FAI在藻华提取过程中的阈值并没有有效的方法。

综上所述，本文针对提取藻华时FAI指数阈值难以确定这一问题，提出浮游藻类指数线性拟合模型(FAI linear fitting model，FAI-L)，先利用坡度计算确定NDVI阈值，再通过建立NDVI与FAI线性回归方程确定FAI阈值的方法，该方法能够较好地获取巢湖水域藻华的空间分布情况，且具有较高的精度和适用性，可为藻华提取光谱敏感指数的阈值确定提供新方法。

1 研究区域与数据源 1.1 研究区域概况

巢湖(31°25′~31°43′N，117°17′~117°51′E)是中国五大淡水湖之一^[23]，同时也是安徽省内最大的淡水湖。湖岸线约为170 km，面积约为780 km²，湖水平均水深3 m，最大水深6.78 m，主要入湖河流为店埠河、南淝河、派河、丰乐河、杭埠河、白天石河、柘皋河等。巢湖蓄含大量淡水资源，是江淮地区重要的生态基础^[24]。近年来随着我国城市化进程不断加快，特别是巢湖临近城市——合肥市的快速发展，大量生活、工业污水排入巢湖，导致湖泊富营养化严重，藻华现象不断发生^[25]，严重影响周边生态系统平衡和居民用水安全。

1.2 数据及预处理

本研究使用Landsat8-OLI和Sentinel-2两种遥感影像，Landsat8共4幅影像，Sentinel-2共6幅影像。其中Landsat8数据通过美国地质勘探局USGS(https://earthexplorer.usgs.gov/)网站下载。Sentinel-2通过欧洲委员会联合欧洲太空总署ESA(https://scihub.copernicus.eu/dhus/#/home)下载。影像产品云量较少，符合研究要求。遥感影像预处理环节包括：辐射定标、大气校正、研究区域裁剪等环节。为了验证研究方法的准确性，选取不同卫星的同时期影像进行对照对比分析，选取影像的类型和时间见附表Ⅰ。

2 研究方法 2.1 藻华提取流程

本文利用遥感数据提取巢湖区域藻华范围(提取流程见附图Ⅰ)，首先通过相关网站下载原始影像数据，其次需要进行相应的预处理，包括辐射定标、大气校正、几何校正等环节，其中Landsat8数据采用ENVI中FLASSH模块进行大气校正，Sentinel-2通过SNAP中Sen2Cor插件进行大气校正，最后通过NDVI和FAI-L提取藻华分布范围，由于巢湖沉水植物主要分布于湖岸附近，本文忽略由于沉水植物对提取结果的干扰。NDVI阈值采用“陡坡效应”原理计算得出。FAI阈值通过浮游藻类指数线性拟合模型(FAI-L)确定，得到提取的藻华范围。并对提取结果进行精度评定分析、方法适用性分析，探究该方法对于不同数据中的藻华提取精度。

2.2 波段选择

由于本文采用NDVI和FAI-L两种方法提取藻华面积，需要使用红光，近红外和短红外进行波段运算。由于Sentinel-2在近红外范围内波谱较宽(共有13个波段覆盖了可见光、近红外、短红外^[26])，为了尽可能降低由近红外波段中心波长偏差导致的藻华提取误差，本文以Landsat8卫星近红外波段为参照(图 1)，通过对比光谱反射率选取合适波段^[27]，最后选取Sentinel-2的B8A(865 nm)波段作为近红外波段参与NDVI与FAI-L计算。

2.3 基于NDVI获取藻华空间分布方法

在对藻华水体识别的众多方法中归一化植被指数(NDVI)是最常见的方法之一^[28]，具有良好的稳定性^[29]。其公式表达为：

$ N D V I=\left(\rho_{\text {nir }}-\rho_{\text {red }}\right) /\left(\rho_{\text {nir }}+\rho_{\text {red }}\right) $

(1)

式中，ρ_nir为近红外波段，分别对应Landsat8中的B5波段、Sentinel-2中的第B8A波段，ρ_red为红光波段分别对应Landsat8中的B4波段、Sentinel-2中的第B4波段。

当水体中藻类浓度升高时，水体在红光波段表现出很强的吸收特性，使NDVI值趋向于正值，而清洁水体则刚好相反，使NDVI值趋向于负值。为了提高提取精度，本文采用坡度与统计分析的方法确定NDVI的阈值^[30]。首先，为防止高藻华区域对坡度分类造成影响，先将NDVI>0.2的值判定为高浓度藻华区域^[31]，再对NDVI < 0.2的区域使用ENVI软件中的Slope工具进行坡度计算，得到NDVI坡度图。其次，在ArcGIS中对像元坡度进行统计分析，按照自然间断点分级法(Jenks)将坡度像元分为两高坡度和低坡度两类^[20]，低坡度像元所对应的为清洁水体范围，高坡度像元所对应的NDVI值为藻华的提取范围。最后，将提取的NDVI均值作为藻华阈值，获得藻华的空间分布情况。

2.4 基于FAI-L获取藻华空间分布方法

FAI指数是由Hu^[17]提出的藻类敏感指数。由于水体对红光波段具有较高的反射率，而在近红外波段和短红外波段具有较强的吸收能力。当水体中藻华浓度上升时，藻华水体对红光波段的吸收能力增强，使藻华水体能够与清洁水体区分开来。因此FAI指数可用于藻华的提取，其公式表达为：

$ F A I=\rho_{\text {nir }}-\rho_{\text {nir }}^{\prime} $

(2)

$ \rho_{\text {nir }}^{\prime}=\rho_{\text {red }}+\left(\rho_{\text {swir }}-\rho_{\text {red }}\right) \frac{\lambda_{\text {nir }}-\lambda_{\text {red }}}{\lambda_{\text {swir }}-\lambda_{\text {red }}} $

(3)

式中，ρ_red、ρ_nir、ρ_swir分别对应红色、近红外、短红外波段，λ_red、λ_nir、λ_swir分别对应红光、近红外、短红外波段的中心波长。目前关于FAI指数的阈值判定往往依据经验模型^[32]，缺乏相应的数据支撑，对提取结果造成一定影响。因此，本文提出FAI-L，假设：若NDVI和FAI都能够较好地区分藻华水体和清洁水体，则同一遥感影像的NDVI和FAI值可能存在一定的线性关系，那么可以通过回归分析法建立两者的关系模型，将NDVI阈值带入线性回归方程中，通过NDVI阈值确定FAI阈值。具体流程为：1)通过坡度计算(slope)得到NDVI的阈值；2)在ArcGIS软件中均匀生成500个随机点(随机点分布如图 2)，将同一点处的NDVI值与FAI值分别赋予随机点；3)剔除数据异常值后，确定NDVI与FAI回归参数，建立回归方程；4)将NDVI阈值带入回归方程中计算得到FAI阈值。

3 结果及分析 3.1 基于NDVI的藻华提取

目前基于NDVI指数的藻华提取已经得到广泛的应用^[33]。本文首先使用Landsat8数据，采用Slope坡度计算的方法确定NDVI指数的藻华提取阈值，通过计算得到20180731、20200805、20201024、20210808的阈值分别为-0.051、-0.019、-0.157、0.004，代入相对应的阈值后获得相对应时相的藻华空间分布(分布结果如图 3)。为了判断提取结果的准确性，本文采用目视解译的方法对提取结果进行精度评定，以近红、红和绿波段合成的假彩色影像作为参考(假彩色影像见附图Ⅱ、附图Ⅲ)，由于藻华提取出现的误差主要存在于藻华像元与水体像元的边界处，因此本文选取的样本主要位于两者的边缘区域，对每幅影像抽取超过1000个样本，将选取样本分为藻华像元和水体像元，然后与提取结果相比较，得到ture water(TW)、false water(FW)、ture cyanobacterial(TC)、false cyanobacterial(FC)，TW表示样本点判断为水体像元且提取结果也显示为水体像元的像元数，FC表示样本点判断为水体像元但提取结果显示为藻华像元的像元数，TC表示样本点判断为藻华像元且提取结果也显示为藻华像元的像元数，FW表示样本点判断为藻华像元但提取结果显示为水体像元的像元数。计算它们的用户精度(P_u)、总体分类精度(P₀)和Kappa系数(K)。对应的计算公式如下：

$ P_0=\frac{T C+T W}{T C+F W+F C+T W} \times 100 \% $

(4)

$ P_1=\frac{(T C+F W)(T P+F C)+(F C+T W)(F W+T W)}{(T C+F W+F C+T W)^2} $

(5)

$ K=\frac{P_0-P_1}{1-P_1} $

(6)

$ P_{\mathrm{u}}=\frac{T C}{T C+F W} $

(7)

式中，P₁为Kappa系数过程变量，K为Kappa系数。

提取结果如表 1所示，通过不同数据源的精度分析可以看出Landsat8提取结果的总体精度平均为91.72 %，为了探究造成误差的原因，通过观察20180731 NDVI提取结果可以发现，Landsat8影像提取结果中东北角存在大面积藻华区域，观察假彩色合成影像发现影像东北角对应的地区为较大云块，认为NDVI提取结果存在将云块阴影误判为藻华，导致提取结果差距较大。同样地，对于20200805 Landsat8中的云块也被误判为藻华区域，可见云层对于NDVI指数有较大影响，容易导致藻华区域的误判。

3.2 基于FAI-L的藻华提取结果

本文提出建立NDVI与FAI指数的线性回归函数关系，建立浮游藻类指数线性拟合模型(FAI-L)，该模型是通过NDVI阈值范围，确定FAI阈值范围的方法。首先对Landsat8数据进行FAI阈值计算。由于各类指数的阈值范围对于不同遥感数据甚至同一数据的不同时间，均存在一定的差异。为了提高阈值的精确度，本文所用影像均建立各自的线性回归函数，得到拟合曲线如图 4所示，将NDVI阈值代入拟合函数中计算得到各个影像FAI阈值。通过拟合计算分别得到20180731、20200805、20201024和220210808 Landsat8影像的NDVI和FAI线性回归方程(各自相关系数R均大于0.8，决定系数R²大于0.65)。拟合效果较好，可以满足实验需求。将各个影像所对应的NDVI阈值带入线性回归方程中，计算得到FAI阈值，通过对应的阈值，得到FAI-L提取藻华的空间分布，结果如图 5所示。

为了验证FAI-L提取藻华空间分布的准确性，本文以假彩色合成图像为参照，通过目视解译选取藻华和水体像元(样本数量>1000)作为样本与提取结果相对比，对各影像数据的提取结果进行精度评定，评定结果如表 2所示。从表 2可以看出，对于Landsat8数据FAI-L在提取藻华的平均精度上可以达到97.16 %，相对于NDVI的提取精度(91.72 %)提升了5.44 %。可见本文采用构建NDVI和FAI回归方程的方法，能有效提高藻华提取的精度。

3.3 FAI-L的适用性分析

为了验证本文中采用的FAI阈值法的适用性，采用Sentinel-2数据进行验证，以20180731、20200725、20200824、2021073、20210819、20210918共6景影像为数据源，提取相关区域的藻华空间分布情况。首先使用NDVI阈值法提取相关区域的藻华分布，同样采用坡度计算的方法计算出NDVI阈值，计算后得到6景Sentinel-2影像的NDVI阈值见表 3。通过阈值获得各影像的藻华分布情况(图 6)。随后通过构建Sentinel-2数据NDVI与FAI线性回归函数如图 7所示(各自相关系数R均大于0.8，决定系数R²大于0.7)，代入NDVI阈值后计算得到20180731、20200725、20200824、2021073、20210819、20210918的FAI阈值见表 3。通过对应阈值得到FAI提取的藻华空间分布，与NVDI提取结果相对比，结果如图 8所示。

对Sentinel-2影像数据FAI-L的提取结果进行同样的精度评定，结果如表 4所示。可以发现该方法对于分辨率更高的Sentinel-2数据有着更高的适用性，对于Sentinel-2数据采用NDVI提取结果的平均精度为93.25 %，采用FAI-L的提取结果平均精度为97.10 %。FAI-L的提取精度较NDVI精度提高了3.85 %。与Landsat8数据的提取结果相比，采用Sentinel-2数据时NDVI提取精度提高了1.35 %，而两者的FAI-L的提取精度相近。因此认为该方法同样也适用于分辨率较高的影像数据。

3.4 NDVI提取结果分析

20180731的Landsat8和Sentinel-2数据为同一天拍摄，但两者的藻华提取结果存在一定差异。推测可能是由于卫星影像拍摄时间不一致，或者是由于卫星分辨率不同而造成的差异。为了探究造成差异的原因，本文查阅Landsat8和Sentinel-2影像的成像时间分别为当天10点42分50秒和10点45分39秒，两者相差的时间为2分49秒，时间相对较短，因此认为由于时间差异造成的误差比较有限。以分辨率最高的Sentinel-2假彩色影像为参考，采用局部放大观察的方法探究NDVI的提取情况，结果如图 9所示。从图 9中d~f、g~i和j~l对比可以看出，对于分辨率较高的Sentinel-2数据在藻华的提取上更加精细，对于较为细小的藻华条带可以完整的提取，如1#、2#和3#三个区域的藻华斑块，其边缘提取情况也与假彩色图像更加吻合。可见分辨率较高的Sentinel-2对细小的藻华有较强的识别能力。分辨率较高遥感影像对细小的藻华识别能力可能更强，提取的藻华范围可能与实际更相符。

3.5 NDVI和FAI-L两类方法的对比分析

通过对比NDVI与FAI-L提取的藻华面积(表 4)发现两者在提取面积上非常接近(相对误差均值在3 % 左右)，说明FAI-L提取藻华的有效性。在提取结果的精度分析中可以看出FAI-L的总体精度平均值为96.15 %，而NDVI阈值法的总体精度仅为91.72 %。可见FAI-L在提取精度上较NDVI提高了约5 %。且无论是Sentinel-2或Landsat8数据，相较于NDVI指数提取精度均有提高，因此本文认为FAI-L在藻华提取精度上优于NDVI阈值法。将20180731 Landsat8和Sentinel-2的NDVI和FAI-L提取结果与相应的假彩色影像进行对比，结果如图 10。从假彩色合成图像可以观察到两者东北角均存在大面积云层(图 10中红色矩形区域)，而NDVI将部分云层误判为藻华水体，而FAI-L并未误判。通过对矩形中藻华面积进行计算得到Landsat8和Sentinel-2误提取面积分别为10.43、18.24 km²，可见由于存在云层误判的情况导致提取结果误差较大。因此FAI-L提取藻华可以有效排除云层的干扰，提高提取精度。

4 结论

本文针对NDVI和FAI提取藻华时存在的阈值难以确定的问题，提出了浮游藻类指数线性拟合模型(FAI-L)。以Landsat8数据为基础，通过建立NDVI与FAI回归方程确定FAI阈值，可以有效提高藻华提取精度。并且分析了FAI-L方法对Sentinel-2数据的适用性及其优点，研究得到的主要结论如下：

1) 对于Landsat8数据，相较于NDVI阈值法，采用FAI-L方法在提取精度上有明显提高，精度由91.72 % 提高至97.16 %，精度提高了5.44 %。

2) FAI-L方法表现出对高分辨率卫星数据如(Sentinel-2)有较强的适用性，对藻华的提取精度可以达到97.10 %。

3) 由于部分影像存在云层的干扰，NDVI阈值法难以将云层和藻华区域区分开，容易造成误判。而本文采用的FAI-L方法可以有效避免云层干扰，当有一定云雾干扰时，采用FAI-L可以有效提高提取精度。

5 附录

附表Ⅰ和附图Ⅰ~Ⅲ见电子版(DOI: 10.18307/2023.0416)。