近岸水体因其海陆交互作用强烈的特殊地理位置，与周边的社会经济构成了多元化的复合生态系统，随着环湾地区社会经济的发展，水体富营养化现象日益严重^[1-3]。Chl a浓度常作为表征水体富营养化程度的重要指标。传统的水质监测是通过采集表层水样、过滤、萃取以及分光光度计分析来测定Chl a浓度。这种方法虽然数据精度相对较高，但耗时耗力，以点带面的不足性明显，不能实现大区域的实时监测。遥感技术以其快速实时、监测范围广、频率高等优势，在各类水体的水质动态监测中发挥了极其重要的作用^[4-6]。

水质遥感反演主要是通过建立特定波段的反射率特征和Chl a浓度之间的关系模型，从而反演Chl a浓度。水体按其光学性质可分为Ⅰ类水体和Ⅱ类水体^[7]。Ⅰ类水体的光学特性主要由浮游植物及其伴生物决定，典型的Ⅰ类水体是大洋开阔水体。Ⅱ类水体的光学特性相对复杂，除了浮游植物，还受到非色素悬浮物、黄色物质等的影响，这类水体主要位于近岸、河口等受陆源物质排放影响较为严重的地方。因此，对Ⅱ类水体而言，遥感反演的应用和精度都存在更大的不确定性^[8]。国外学者最早开拓了以经验方法为主的水质定量遥感反演，研究表明，在可见光和近红外波段组合的反射率与Chl a浓度有较好的拟合关系^[9-10]。国内学者马荣华^[11]等的研究结果表明Chl a浓度与ETM第三波段与第一波段的比值组合有较高相关性。90年代后国外学者开启了以半经验半分析方法为主的Chl a反演算法。Tedesco L P^[12]发现在水质复杂的近岸水域，三波段的遗传算法能较好的反演Chl a浓度，决定系数R²在0.91以上，均方根误差RMSE (root mean square error)小于5 mg/m³。Giorgio D O^[13-15]等人基于生物光学模型表明MERIS(medium resolution imaging spectrometer)中心波长为708 nm、665 nm、778 nm的波段反射率可以用于陆海交互作用强烈的水体中Chl a的反演，并获得了较好的效果。Stefani N^[16]等结合MERIS数据，在比斯开湾建立了反演Chl a的波段比值算法，R_rs(remote-sensing reflectance)(510)/R_rs(560)，决定系数R²达到0.681。国内学者周冠华^[17]和郝艳玲^[18]通过最佳波段组合反演Chl a浓度，决定系数R²分别达到0.8358和0.9。韩秀珍^[19]和张鸽^[20]基于FY-3A/MERSI(medium resolution spectral imager)分别建立了太湖水体Chl a和悬浮物浓度反演模型，决定系数分别为0.60、0.72和0.89、0.78。

胶州湾水质为典型的Ⅱ类水体，Chl a浓度的时空动态规律和发展趋势存在较大不确定性，不利于水质的动态监测和治理。因此，本研究以Landsat-8、FY-3A和HJ-1B遥感影像为数据源，结合实测的Chl a浓度，在分析水体光谱特征的和综合前人研究基础上，建立能较好适应于胶州湾Chl a浓度反演的半分析模型，同时，对其他的Ⅱ类水体的相似研究有一定的借鉴意义。

1   材料与方法

1.1   研究区域

胶州湾(35°58′N~36°18′N，120°04′E~120°23′E)位于山东半岛的南岸，地处亚洲大陆边缘，是典型的温带半封闭型海湾(见图 1)，面积约为423 km² ，平均水深约为7 m，湾口及湾中部的水深达30 m以上。目前胶州湾入海河流有十几条，如大沽河、洋河、海泊河、李村河、娄山河等，其中娄山河、李村河和海泊河已成为几个主要的排污口。生态环境受自然和人类活动的双重影响，自然因素主要有东亚季风、黄海水团系统等，而人类活动主要包括陆源排污及围海养殖等。

图  1  研究区域和实测采样点分布

Fig.  1  Map of the study area and in situ sampling sites

下载: 全尺寸图片 幻灯片

1.2   采样点分布及水质分析

本研究分别于2016年2月9日、5月2日、8月3日、11月23日在胶州湾海域进行水样采集，严格按照《海洋监测规范》规定执行，采样点分布如图 1所示。实验室Chl a浓度的测定通过过滤、萃取、分析和测量，采用分光光度法(GB 17378.7-2007)。Chl a浓度年际变化呈现双峰型，在2月和8月浓度值出现高峰，见表 1。

表  1  Chl a浓度统计结果

Tab.  1  Statistical results of Chl a concentration

下载: 导出CSV 

| 显示表格

1.3   实测光谱数据

水体光谱数据采用的是水面测量法。使用美国ASD(analytical spectral devices)公司生产的FieldSpec-FR野外光谱仪，其扫描范围是325~1075 nm，采样间隔为1 nm，仪器的视场角为25°。采样时间为9:30~15:40，采样当天微风。为避免太阳的直射以及船体的阴影会对测量产生影响，测量时背对太阳，每个采样点测量10条光谱数据，以减少光照和仪器稳定性对测量的影响。待测物理量：总辐射亮度L_sw、天空漫散射光亮度L_sky、标准灰板辐射亮度L_p。水体遥感反射率R_rs按照下式进行计算：

式中：r为水-气界面反射率，其影响因素有风速、太阳高度角和观测角度等，取值为0.026~0.028；ρ_p为标准板反射率，测量中选用的是30%的标准板；π取3.1415926。

为了消除不同时间点、不同站位、不同站位光谱测量时外界环境因素变化产生的影响，故对水体野外实测光谱曲线数据进行归一化处理。水体遥感反射率数据归一化处理的公式为:

式中：R_N(λ_i)是波长λ_i处光谱曲线归一化值；R(λ_i)野外实测光谱曲线的反射率；n是400~900 nm波段之间的波段数。

1.4   卫星数据预处理

根据遥感影像的同步性要求，本研究选取2016年的卫星影像过境时间分别为：Landsat-8/ OLI(operational land imager，陆地成像仪)，2月9日、4月29日、8月3日和11月23日；FY-3A/MERSI(中等分辨率光谱成像仪)，2月9日、5月3日、8月4日和11月23日；HJ-1B/CCD(charge-coupled device可见光相机)，2月9日、5月1日、8月3日和11月24日。对所获取影像均进行预处理工作。首先进行控制点精校正，采取二次多项式，使误差控制在0.5个像元之内；其次进行辐射定标，将传感器记录的数字量化值DN(digital number)转换成绝对辐射亮度值(辐射率)，以精确反演地物特征；最后进行大气校正，目的是消除大气和光照等因素对地物反射的影响，以反演地物真实反射率，采用6S模型进行较正。

2   结果与讨论

2.1   水体光谱曲线归一化处理

由公式(1)(2)得出2016年4个月份的水体归一化光谱曲线，剔除异常光谱，每个采样点的光谱为10次光谱的平均。图 2为5月份归一化前后水体的光谱曲线图。

图  2  归一化前后光谱曲线

Fig.  2  Spectral curves and normalized processing

下载: 全尺寸图片 幻灯片

胶州湾水体的反射率呈现出内陆水体光谱特征，光谱带中Chl a的光谱特征明显，即在420 nm~500 nm、675 nm有Chl a的吸收峰，在700 nm左右有Chl a的反射峰。波长在400 nm~500 nm范围内反射率较低，是由于藻类色素和黄色物质在该范围内的强烈吸收作用；波长在500 nm~590 nm范围内的反射峰是藻类光合色素的弱吸收和悬浮颗粒物的散射作用共同形成的，该范围内2月份的光谱曲线变化趋势较其余月份波动性较大；波长在670 nm附近为Chl a的吸收峰，反射率出现谷值；波长在685 nm~715 nm范围内反射峰的出现是含藻类水体最显著的光谱特征，可作为判定水体是否含有藻类Chl a的依据。

2.2   敏感性波段选取

为了更好地构建胶州湾的反演模型，提高反演精度，对实测Chl a浓度和同步水体光谱曲线之间进行相关性分析，选取敏感波段。400 nm~900 nm波长范围的归一化实测遥感反射率与Chl a浓度的相关系数如图 3所示。可知，2月份两个较高正相关波段(704 nm和657 nm)及较高负相关波段(400 nm~475 nm)能较好预测此时的水体Chl a浓度；5月份，在588 nm附近，水体反射率与叶绿素浓度a呈现较高负相关，波长大于700 nm出现了较高的正相关；8月份，波长在650 nm和700 nm附近有较高的正相关性，呈现较好负相关性的波段位于475 nm~510 nm附近；11月份，波长在410 nm~480 nm附近水体遥感反射率与Chl a浓度呈现良好的正相关，在700 nm附近的负相关系数最高。由此可以通过选取较高正相关波段和较高负相关波段对胶州湾水体Chl a浓度进行反演。

图  3  Chl a浓度与归一化遥感反射率的相关系数

Fig.  3  Correlation coefficient of chlorophyll-a concentration and normalized R_rs

下载: 全尺寸图片 幻灯片

2.3   基于HJ-1B反演模型的建立与分析

HJ-1B/CCD的波段设置：Band1(Blue)，0.43 μm~0.52 μm；Band2(Green)，0.52~0.60 μm；Band3(Red)，0.63 μm~0.69 μm；Band4(NIR)，0.76 μm~0.90 μm。根据2.2分析可知，为了较好地模拟Chl a浓度值，2月份选择反演模型CCD Band 3/CCD Band 1，5月份选择反演模型CCD Band 4/CCD Band 2，8月份选择反演模型CCD Band 3/CCD Band 1，11月份选择反演模型CCD Band 1/CCD Band 4。从2016年的野外实测数据中随机选取18个点建立模型，预留8个点用来检验模型的精度，预留验证点(见图 1)均匀覆盖研究区域，具有验证模型的可行性。同时对每一个模型采用五种函数：线性、指数、幂函数、对数和一元二次多项式算法进行拟合，结果如表 2。

表  2  基于HJ-1B的反演模型和精度检验

Tab.  2  Inversion model based on HJ-1B and accuracy test

下载: 导出CSV 

| 显示表格

据表 2可知，基于HJ-1B/CCD Chl a浓度反演的模型中，一元二次多项式算法在4个月份中的决定系数均最高，除8月份的线性函数和5月份的幂函数达到0.6以外，其余函数的拟合效果较差。2月份最高决定系数R²=0.63，RMSE=1.22 mg/m³；5月份的最高决定系数R²=0.67，RMSE=1.67 mg/m³；8月份的最高决定系数及对应的均方根误差结果都较好，R²=0.69，RMSE= 0.90 mg/m³；11月份的最高决定系数R²=0.66，RMSE=1.10 mg/m³。综上所述，该生物光学模型反演的相关系数整体精度值偏低，其主要用来监测水体的蓝色波段波长范围(0.43 μm~0.52 μm)较Landsat系列相对较宽，影响了水体的反演精度，不能较好的适用于胶州湾Chl a浓度的反演研究，但对其他的相似研究有一定的参考价值。

2.4   基于FY-3A的反演模型的建立与分析

FY-3A/MERSI光谱分辨率高，共20个通道，其中6-14通道可用于水色遥感。从可见光到近红外光谱范围内的四个波段中，其中心波长为470 nm，550 nm，650 nm的波段靠近Chl a的敏感波段。目前FY-3A/MERSI数据在内陆水体的水质参数反演方面的研究不多，仍处于探索和发展。国际上常用的MERIS半分析反演模型，为Dall’Olmo^[21]等提出的三波段模型：(=671 nm，=710 nm，=740 nm)。本研究基于对胶州湾水体敏感波段的分析和MERSI的波段设置，构建了胶州湾水体的波段比值模型。2月份，=679 nm，=682 nm，=704 nm；5月份，=671m，=715 nm，=754 nm；8月份，=673 nm，=704 nm，=715 nm；11月份，=682 nm，=704 nm，=720 nm；结果如表 3。

表  3  基于MERSI的反演模型和精度检验

Tab.  3  Inversion model based on MERSI and accuracy test

下载: 导出CSV 

| 显示表格

据表 3可知，MERSI的反演结果中，一元二次多项式算法的决定系数最高，R²为0.69。除5月份外，相同月份的三波段一元二次多项式拟合的决定系数均高于两波段比值模型。2月份的最高决定系数R²=0.65，RMSE=1.24 mg/m³；5月份的最高决定系数R²=0.59，RMSE=1.53 mg/m³；8月份的最高决定系数R²=0.69，RMSE=1.40 mg/m³；11月份的最高决定系数R²=0.60，RMSE=0.87mg/m³。鉴于前人研究^{[13-16, 19-20]}的结果，该模型反演的整体决定系数偏低，不能适用于胶州湾叶绿素浓度的反演研究。决定系数出现较低的原因可能是胶州湾水体受海陆交互作用影响较大，随着环湾产业区的发展，污染日益加重，导致水质较复杂。同时，MERSI的空间分辨率较大，降低了水质复杂的Ⅱ类水体的反演精度。

2.5   基于Landsat-8反演模型的建立与分析

将目前主流的半分析模型应用于Landsat-8/ OLI卫星数据^[22-24]，同时结合2016年胶州湾水体的光谱特性，选取前五个波段：Band1(Coastal)，0.43 μm~0.45 μm；Band2(Blue)，0.45 μm~0.51 μm；Band3(Green)，0.53 μm~0.59 μm；Band4(Red)，0.64 μm~0.67 μm；Band5(NIR)，0.85 μm~0.88 μm，据此调整了二波段模型、三波段模型。2月份反演模型为Band4/Band2、[(Band2)^-1-(Band3)^-1]* Band4；5月份反演模型为Band5/Band3、[(Band3)^-1-(Band4)^-1]* Band5；8月份反演模型为Band4/Band3、[(Band2)^-1-(Band3)^-1]* Band4；11月份反演模型为Band1/Band4、[(Band4)^-1-(Band2)^-1]* Band1。同时对每一种模型采用五种函数：线性、指数、幂函数、对数和一元二次多项式算法进行拟合。结果如表 4。

表  4  基于Landsat 8的反演模型和精度检验

Tab.  4  Inversion model based on Landsat 8 and accuracy test

下载: 导出CSV 

| 显示表格

据表 3可知，基于Landsat半分析模型反演的整体效果较好，五种函数类型中，除对数函数外，其余函数类型中90.625%的决定系数R²均在0.6以上，能适应于胶州湾叶绿素浓度的反演。与3.1和3.2的结果相似，一元二次多项式算法在两种模型中拟合的效果均最好。除5月份最高决定系数R²=0.77外，其余最高决定系数R²均在0.8以上。2月份两波段模型反演优于三波段模型，5、8和11月份均是三波段模型优于两波段模型，与2月份的光谱特征有关，在500 nm~590 nm范围内的光谱曲线波动性较大，因此红光波段(Band3:0.53 μm~0.59 μm)不适用于2月份Chl a浓度的反演。2月份最高决定系数R²=0.83，RMSE=1.07 mg/m³；5月份最高决定系数R²=0.80，RMSE=1.0 mg/m³；8月份最高决定系数R²=0.78，RMSE=0.89 mg/m³；11月份最高决定系数R²=0.86，RMSE=1.13 mg/m³；故能较好的适应于胶州湾叶绿素浓度的反演研究。图 4为4个月份Chl a实测浓度值与最优模型反演结果图。为了检验上述最优模型的科学合理性，对模型进行误差敏感性分析，将预留的8个实测点浓度值与对应算法预测浓度值进行对比分析，以检验模型算法的可靠性，结果如图 5所示。

图  4  Chl a浓度与波段反射率相关关系

Fig.  4  Correlation between chlorophyll-a concentration and band reflectance

下载: 全尺寸图片 幻灯片

图  5  Chl a实测浓度与模型预测值相关关系

Fig.  5  Correlation between in-situ concentration of chlorophyll-a and model predictive value

下载: 全尺寸图片 幻灯片

2.6   讨论

胶州湾水体2月份Chl a的较高正相关波段在704 nm和657 nm附近，较高负相关波段位于400 nm~475 nm之间，5月份波长大于700 nm时呈现较高的正相关性，波长在588 nm附近呈现较高的负相关性；8月份较高正相关波段在650 nm和700 nm附近，较高负相关波段位于475~510 nm之间；11月份波长在410~480 nm附近呈现良好的正相关性，波长在700 nm附近的负相关性最高。据此调整的HJ-1B/CCD反演模型为：2月份CCD Band 3/CCD Band 2，5月份CCD Band 4/CCD Band 2，8月份CCD Band 3/CCD Band 1，11月份CCD Band 1/CCD Band 4，决定系数R² < 0.7，整体偏低，不能较好的适用于胶州湾Chl a浓度的反演研究，但对其他的相似研究有一定的参考价值；FY-3A/MERSI反演模型为，决定系数R² < 0.69，整体偏低，不能适用于胶州湾复杂水质中Chl a浓度的反演研究。Landsat-8/OLI最佳波段组合为：2月份Band4/Band2，5月份[(Band3)^-1-(Band4)^-1]* Band5，8月份[(Band2)^-1-(Band3)^-1]* Band4，11月份[(Band4)^-1-(Band2)^-1]* Band1，决定系数R²>0.80，预测值与实测值的决定系数R²>0.86，能较好的适用于胶州湾Chl a浓度的反演研究。

3   结论

(1) 本研究以Landsat-8/OLI、FY-3A/MERSI和HJ-1B/CCD卫星影像为数据源，并结合2016年实测的数据，对胶州湾Chl a浓度进行遥感反演。通过分析实测的光谱数据，筛选出对应卫星适合反演Chl a浓度的波段组合，建立了能较好的适应于胶州湾复杂水体中Chl a浓度的遥感反演模型，为研究胶州湾Chl a浓度的时空变化规律和水质监测奠定了基础。

(2) 由于本研究采用的实测样本数据相对较少，不能明确分析出反演结果与卫星影像的空间分辨率和光谱分辨率的之间的关系以及各自的影响权重，模型的代表性及稳定性还需要大量的实测数据进一步验证。在以后的工作中可以加大采样点数量，使各模型中趋于更加稳定，即可以在没有同步地面数据的情况下通过卫星遥感图像和模型来估测当时的Chl a浓度。