基于GF-1/WFV的钱塘江叶绿素a和总悬浮物浓度遥感估算

， ， ，，

(1.浙江水利水电学院 测绘与市政工程学院，杭州 310018；2.浙江工商大学 旅游与城乡规划学院，杭州 310018；3.浙江省水文局，杭州 310024； 4.浙江水利水电学院 水利与环境工程学院，杭州 310028)

1 研究背景

水体中的叶绿素a和悬浮物是重要的水质参数。叶绿素 a 是浮游植物或藻类植物中最丰富的色素，其含量大小在一定程度上反映了水体富营养化的程度[1]。悬浮物是指悬浮在水中的固体物质，包括不溶于水的无机物、有机物、泥沙、黏土和微生物等，是我国近岸海水环境和湖泊水环境质量检测的重要参数之一[2-3]。常规的水质参数人工监测方法费时、费力，且由于河流水体流动性大，采样点的数目和次数有限，监测结果难以代表整个水体中水质指标的分布状况。利用遥感方法监测水体水质范围广、速度快、成本低，可以反映水质参数在空间和时间上的分布情况和变化，便于长期动态监测。

目前国内外主要是利用基于中低分辨率遥感数据如MODIS，MERIS，Landsat/TM，HJ1 /CCD和高分辨率遥感数据如SPOT，GF-1等对河流的典型水色参数叶绿素 a浓度和悬浮物浓度进行反演模型研究[4-6]。然而，这些研究主要关注面积较大的内陆湖泊或海洋等水体，对于宽度较小、水体流动性大的河流水体关注较少；另外针对与城镇人口生活关系密切的河流水体，必须使用高空间分辨率的遥感数据进行水质参数估算，才能快速监测河流水质参数的空间分布变化特征。2013 年 4 月，我国自主研制的第一颗民用高分辨率卫星高分一号( GF-1)成功发射，GF-1卫星搭载了2台2 m分辨率全色/8 m分辨率多光谱相机和4台16 m分辨率多光谱相机，光谱信噪比高，为宽度较小、水体流动性大的河流水体进行典型水色参数遥感估算提供了重要的硬件基础。

钱塘江水系是浙江省最大的水系，是浙江省重要命脉，钱塘江水环境污染问题会严重制约经济社会的可持续发展[7]。基于高分一号(GF-1)16 m分辨率的GF-1/WFV传感器数据和地面实测数据，利用地面实测光谱反射率数据模拟的WFV波段组合，构建水体叶绿素a和总悬浮物浓度遥感模型，对钱塘江流域杭州段水体水质进行遥感估算，以期探索利用高分一号卫星数据遥感监测河流水色参数，也为后续基于高分遥感影像的河流水质监测的理论和方法研究提供重要的参考依据。

2 研究区概况和数据源处理

2.1 研究区与样点

研究区域为钱塘江杭州段之江水文站断面，采样区域长2 km，宽1 km，各样点之间间隔300～500 m，低潮位时水深约3 m，测量样点15个，样点分布如图1所示。测量时间为2016年12月7日9：00—16：00，间隔0.5 h取样测量一次,测量内容包括光谱测量和水样采集。

图1 研究区及采样点分布Fig.1 Distribution of sampling points in the study area

2.2 实验数据获取与处理

野外光谱测量使用美国分析光谱仪器公司的ASD(Analytical Spectral Devices)野外光谱辐射仪(ASD FieldSpec.)。为减少水体镜面反射和船体自身阴影的影响，测量采用二类水体水面以上光谱测量方法[8]，测量时距离水面0.5～1.0 m。测量的数据有标准灰板(反射率为30%)、水体和天空光的光谱辐亮度信息，每个对象采集20条光谱数据，选取每个样点重复测量的遥感反射率的中值作为该点的反射率。在测量得到的光谱曲线中，选择350～900 nm的数据计算遥感反射率[8]。

野外采集的水样避光保存，并于当天送到国家海洋局第二海洋研究所卫星海洋环境动力学国家重点实验室进行实验分析，测定叶绿素a和总悬浮物浓度。室内叶绿素a浓度(Chla)使用直径25 mm的Whatman GF/F滤膜进行过滤，丙酮萃取后，使用荧光光度法进行测量。总悬浮物浓度(TSM)使用事先称重的直径47 mm的Whatman GF/F滤膜过滤，并采用常规的干燥、烘烧、称重法进行测定。所有样品均设置2个平行样，2次测量结果的均值作为最终浓度。

由于GF-1/WFV卫星影像的获取时间与实测数据时间存在差异，会影响数据的代表性及模型精度，因此使用实测光谱反射率拟合GF-1/WFV各波段光谱数据进行模型构建。高分一号卫星影像的波段设置为4段：蓝(波长450 ～ 520 nm，)、绿( 520～590 nm)、红(630～690 nm)、近红外(770～890 nm)，其光谱响应函数可以在中国资源卫星应用中心网站上进行获取。根据式(1)，由实测高光谱遥感反射率(波长间隔为1 nm)计算传感器通道i的反射率，传感器光谱通道i的反射率ri为

(1)

式中：λsi为通道i的起始波长；λei为通道i的终止波长；r(λ)为波长λ处的反射率值；φi(λ)为波段i在波长λ处的光谱响应因子。

3 监测数据的反演模型构建与影像处理

3.1 反演模型的构建与评价

基于拟合得到的GF-1/WFV各通道光谱数据与水质参数浓度，建立Chla和TSM反演模型，再将适用性较好的模型应用于GF-1/WFV遥感影像上。实测数据共有15组，选择10组进行反演模型的构建，5组数据进行模型验证。

3.1.1 反演模型构建

基于遥感的水质监测主要是通过研究水体光谱特征与水质参数间的关系构建反演模型。常规的水质参数遥感反演方法主要有经验方法、半经验/半分析方法、分析方法[9]和机器学习算法等[10-11]。其中，利用高光谱或多光谱数据选取特征波段或波段组合，与同步实地测量的水质参数浓度建立半经验/半分析算法，能够较好地反映水体水质参数的光学响应特征，物理意义明确且精度较高，是水体水质参数遥感估算的主要方法，在实际应用中最为广泛。常用传感器的波段设置对比如图2所示,图中括号内的数据为传感器的分辨率。

图2 常用传感器的波段设置对比Fig.2 Comparison of band settings of remote sensing sensors

由于GF-1/WFV传感器的波段设置与Landsat/ETM， HJ1/CCD， SPOT波段设置相似(图2)，因此，借鉴这些传感器的波段组合与模型构建方法，进行叶绿素a浓度(Chla)和总悬浮物浓度(TSM)的反演研究。本文所使用的模型构建方法包括：

(1)基于实测高光谱数据拟合得到的WFV波段，使用穷举法列举各种模型波段组合，包括单波段、基于不同波段两两组合的比值和差值组合等，并针对这些波段组合变量分别进行线性回归、对数线性回归、非线性回归(幂、指数、多项式)拟合建模。

(2)基于单个、2个、3个和4个光谱波段分别进行多元回归拟合建模。

(3)验证已有的特殊波段组合，包括(B1-B4)/(B3-B4)，Ln(B3 + B4)/(B1 + B2)，(B2/B1) + B3，(B3+B4)/B2，(B2×B3)/B1，(B1-B3)/B2，(B4-B3)/(B4+B3)，并根据Chla和TSM的光谱响应特征，参考已有的波段组合构建方式，构建9种新的模型组合，包括(B4-B3)/(B1-B3)，(B4+B3)/(B1+B3)，(B4-B3)/(B1+B3)，(B4+B3)/(B1-B3)，B2/B3+B1，B1/B3+B2，(B2×B1)/B3，(B2+B1)/B3，B1+B2+B3+B4。

最后选择与Chla和TSM相关性最大的多光谱波段组合模型，讨论其模型精度和验证应用精度，经过对比选择适合GF-1/WFV传感器的最优Chla和TSM估算模型。

3.1.2 模型评价

本文使用回归分析构建各波段组合与水质参数浓度间的模型表达式，因此，使用决定系数R2、F检验值和显著性水平p值来评价回归模型的精度和稳健性。基于模型预测值和实测值，使用均方根误差(Root Mean Square Error, RMSE)、平均绝对误差(Mean Absolute Error，MAE)和平均相对误差(Mean Relative Error, MRE)来评价模型的精度，其计算方法如下：

(2)

(3)

(4)

式中:y是测量得到的水质参数浓度;y′是估算出的水质参数浓度;n是样点数目。

3.2 GF-1/WFV影像处理

GF-1/WFV高分影像数据下载自遥感集市网站(http:∥www.rscloudmart.com/)。由于与野外实验同步的WFV高分影像存在30%以上的云层覆盖，影像质量较差，选择与野外实验日期相近的2016年11月28日(农历十月二十九)进行钱塘江杭州段水体叶绿素a浓度和总悬浮物浓度的遥感估算。该期数据的获取时间处于钱塘江低潮位，数据共包括2景WFV2和1景WFV1数据，对其分别进行辐射定标、大气校正、几何纠正和影像拼接等处理，得到钱塘江流域杭州段的WFV数据。

2016年WFV数据的辐射定标公式为

Le(λe)= Gain·DN+ Offset 。

(5)

式中：Le(λe)为转换后辐亮度；Gain 为波段增益；DN 为卫星载荷观测值；Offset为偏置。

数据具体处理步骤如下：

(1)辐射定标。根据绝对定标系数和辐射定标公式，利用ENVI 5.1软件对GF-1/WFV图像进行定标，将卫星载荷观测值转化成辐亮度。GF-1星的WFV相机辐射定标参数，包括4 个波段的增益、偏置参数及定标公式，均可从中国资源卫星应用中心网站上获取。

(2)大气校正。对定标后的大气顶部反射率进行大气校正，使用ENVI 5.1软件的FLAASH模块进行大气校正。将定标后的影像转化成BIL格式，使用10的比例系数将输入光谱数据的单位转换成 μW/ (cm2·sr·nm)，设定传感器高度为645 km，像元大小为16 m，大气模型为Mid-latitude Summer，气溶胶模型为Rural，并根据头文件查找影像获取时间等参数。

(3)几何纠正。使用ERDAS 9.2软件中的Geo-referencing 模块进行几何纠正，以研究区经过几何精纠正的 TM 图像为基准，通过在研究区选取控制点完成对 GF-1/ WFV1和WFV2图像的几何精纠正，纠正后的误差控制在1个像元范围以内。

(4)影像拼接。对预处理后的遥感影像进行拼接，然后使用浙江省杭州市的行政区矢量文件进行掩膜，得到钱塘江杭州段的GF-1/WFV遥感影像图。参考中国科学院计算机网络信息中心地理空间数据云平台上提供的免费Landsat 中国内陆水体信息产品，进行水体范围提取。

4 结果与讨论

4.1 实测数据分析

钱塘江杭州段采集的野外实测光谱和拟合得到的WFV各波段光谱数据(图3)表明，实测光谱反射率呈现出典型的二类水体特征[11]：蓝光波段(350～500 nm)由于水体吸收反射率相对较低，560 nm附近的绿光反射峰和710 nm附近的荧光峰是含藻水体的典型特征，710 nm 附近和810 nm 附近的2个峰值是含悬浮泥沙水体的典型特征。

图3 钱塘江杭州段实测光谱和拟合的GF-1/WFV波段光谱Fig.3 In situ spectral reflectance and fitted GF-1/WFV channel spectra in Hangzhou segment of Qiantang River

图4 钱塘江杭州段水体的叶绿素a和总悬浮物浓度Fig.4 Chlorophyll-a and total suspended matter concentration in Hangzhou segment of Qiantang River

水样的室内分析测定结果表明(图4)，由于采样时间处于低潮位，水质较为清澈，水体中叶绿素a和总悬浮物浓度不高，并且由于样点集中在一个断面附近分布，水质参数浓度变化不大。水体中叶绿素a浓度范围为1.36～7.51 μg/L， 与已有的研究结果较为一致[12]。较多研究认为钱塘江水体的光学特性主要由悬浮物主导，而实测数据表明低潮位时悬浮物的扰动效应较小，叶绿素a与总悬浮物浓度相当。当潮流和径流速度较小时，叶绿素a也会在一定程度上主导河流水体的光学特性，例如2014年钱塘江暴发了严重的蓝藻水华现象。内河段水体中总悬浮物浓度范围为2.43～12.97 mg/L，显著低于钱塘江河口水体的悬浮物浓度[12]，这是由于内河段悬沙含量低于河口，且低潮位采样时风浪对泥沙的扰动作用较小。

4.2 水质参数浓度的遥感估算

4.2.1 叶绿素a浓度遥感估算

基于拟合得到的GF-1/WFV4个波段的光谱数据和实测叶绿素a浓度数据，使用3.1.1节中的波段组合构建叶绿素a浓度反演模型。基于穷举法构建波段组合的模型结果表明，单波段模型中，基于B3波段的多项式拟合结果最优(图5(a))；比值模型中，基于B1/B3组合的多项式拟合结果最优(图5(b))；差值模型中，基于B1-B3组合的多项式拟合结果最优(R2=0.83，RMSE=0.74 μg/L)(图5(c))，并且B1-B3的差值组合拟合结果优于单波段和波段比值组合。

图5 常用的叶绿素a浓度估算模型Fig.5 Commonly used models for estimating chlorophyll-a concentration

此外，基于多个波段数据的多元回归拟合模型也取得了较高的精度，并且使用的波段数量越多，模型精度越高。例如，基于B1,B2,B3,B4这4个波段的多元回归模型精度为R2=0.81，RMSE=0.79 μg/L，MAE=0.70 μg/L，MRE=22.73%，显著高于基于2个波段的多元回归模型。同时，对3.1.1节中已有波段组合和构建的9种新的波段组合进行拟合验证结果表明，特殊波段组合的模型结果与已有的比值或差值模型结果相当，丰富的波段信息和组合变换形式并没有改进模型的精度。上述基于多元线性回归和特殊波段组合的模型结果均低于基于B1-B3的差值组合模型。因此，选择精度最高的差值模型B1-B3作为最终的叶绿素a(Chla)反演模型，其模型形式为

Chla=-367 852×(B1-B3)×(B1-B3)-

467.48×(B1-B3)+ 5.700 6 。

(6)

式中：Chla为叶绿素a浓度； B1为B1波段反射率；B3为B3波段反射率。

式(6)中模型的精度为R2=0.83，RMSE=0.74 μg/L，MAE=0.57 μg/L，MRE=19.67%。使用验证数据集对该模型进行直接验证，结果为RMSE=1.98 μg/L，MAE=1.71 μg/L，MRE=24.97%，表明基于B1-B3差值的二次项模型可以用于钱塘江杭州段的Chla的遥感估算研究。已有研究[13-15]也表明B1和B3波段的比值或差值模型可以用于二类水体Chla的遥感估算，而使用更为广泛的B4/B3波段组合[16-19]在本文的数据中并不适用，可能是由于低潮位时悬浮物浓度对叶绿素a的光谱响应干扰所造成的。

4.2.2 总悬浮物浓度遥感估算

基于拟合得到的GF-1/WFV4个波段的光谱数据和实测总悬浮物浓度数据，使用3.1.1中的波段组合构建TSM反演模型。使用穷举法构建的波段组合的反演模型结果表明，单波段模型中，基于B3波段的多项式拟合结果最优(图6(a))；比值模型中，基于B2/B3组合的线性拟合结果最优(图6(b))；差值模型中，基于B3-B4组合的多项式拟合结果最优(图6(c))。其中，基于B3单波段的二次多项式模型结果最优(R2=0.92， RMSE=1 mg/L)。

图6 常用的总悬浮物浓度估算模型Fig.6 Commonly used models for estimating total suspended matter concentration

对3.1.1节中已有的波段组合和新构建的9种波段组合进行验证，结果表明各波段组合的模型结果与已有的比值和差值模型结果相当，丰富的变换形式对模型精度的改进作用不明显，其中拟合效果最佳的 (B2×B3)/B1，B2/B3+B1和(B2+B1)/B3 组合，拟合精度R2分别为0.88,0.88,0.88，仍低于或约等于常规的比值或差值反演模型的结果。

对比多种模型构建结果，包括基于单波段、波段比值和差值的线性和非线性拟合，以及基于多个光谱波段的多元回归拟合等，发现基于多个波段的多元回归模型精度最高，并且总体呈现出波段数量越多，多元回归模型精度越高的趋势。其中，基于B1，B2，B3，B4这4个波段的多元回归模型为

TSM=-4.516 8+1 926.959×B1-2 420.11×B2+

2 545.432×B3-1 440.59×B4 。

(7)

式中： TSM为总悬浮物浓度； B2为B2波段反射率；B4为B4波段反射率。

式(7)具有最高的反演精度(R2=0.98，F=61.83，p<0.05; RMSE=0.49 mg/L，MAE=0.42 mg/L，MRE=8.72%)。多元回归的反演效果优于基于代数运算所构建的各种模型组合，是一种有效利用传感器多波段信息的方法，然而，多元回归所使用的全部波段信息存在高度相关的现象，如B1与B2相关系数为0.98，B1和B3相关系数为0.89。多元线性回归分析中这种变量的多重相关性常会严重影响参数估计，扩大模型误差，破坏模型的稳健性，从而导致整体的拟合度很大，但个体参数估计值的t统计量却很小，并且无法通过检验[20]。因此，多元回归模型并不适用于本文数据的模型构建。

因此，选择基于B3单波段的多项式模型作为总悬浮物浓度的反演模型，其模型形式为

TSM=57 456×B3×B3-1 041×B3+6.214 1 。(8)

式(8)中模型的精度为R2=0.92，RMSE=1 mg/L，MAE=0.85 mg/L，MRE=17.22%。使用验证数据集对该模型进行直接验证，其结果为RMSE=3.53 mg/L，MAE=2.25 mg/L，MRE=32.95%。针对悬浮泥沙含量高的河口区域，单独使用B4波段进行TSM估算的研究较多[21]，但B4波段在本数据集中应用效果较差(R2<0.6)，可能是由于低潮位的内河段水体中悬浮物受风浪扰动、潮流和径流影响较小，悬浮物浓度较低，波长较短的B3波段更适合对其进行反演。

4.3 基于GF-1/WFV数据的水质参数分布

将前面建立的最佳叶绿素a浓度遥感估算模型(式(6))和总悬浮物浓度遥感估算模型(式(8))代入到经过大气校正的2016年11月28日的GF-1/WFV遥感影像上,进行钱塘江杭州段叶绿素a浓度和总悬浮物浓度估算，得到叶绿素a浓度的空间分布如图7(a)所示，总悬浮物浓度的空间分布如图7(b)所示。

图7 基于GF-1/WFV数据估算的叶绿素a和总悬浮物浓度空间分布Fig.7 Distribution of chlorophyll-a and total suspended matter concentration estimated by GF-1/WFV data

由图7(a)及图1看出，由于低潮位时叶绿素a浓度受潮汐和悬浮泥沙的扰动作用较小，呈现出较高的值，该结果与同期水质监测资料的结果基本一致[22]。从空间分布差异上来看，千岛湖水质较好，叶绿素a浓度最低，钱塘江干流杭州段的叶绿素a 浓度值域范围相差不大，空间差异并不显著。已有研究结果表明，钱塘江叶绿素a 浓度时间差异显著，呈现夏秋季节高、冬春季节低的规律，并且叶绿素a 浓度与温度呈显著正相关，叶绿素a 与总氮量(TN)、总磷量(TP) 之间的相关关系在不同江段有所差异[22]。对比浙江省地表水水质自动监测平台(http:∥wms.zjemc.org.cn/wms/wmsflex/index.html)上提供的监测数据，本研究的叶绿素a分布变化规律与已有的水质调查监测结果中的TN和TP监测结果趋势基本一致，较好地说明了叶绿素a与该水质指标间存在一定的相关性，表明叶绿素a也是反映河流水体富营养化的重要指标之一。

由图7(b)及图1看出，整个钱塘江流域水体总悬浮物浓度比较高，且在空间分布上具有一定差异，其中河口处总悬浮物浓度最高，千岛湖最低，总悬浮物浓度范围是0.23～1 204.5 mg/L，与已有的悬浮物监测结果基本一致[5]。造成钱塘江流域悬浮物浓度空间分布差异的主要原因是人为活动引起的土地利用变化，河口区域聚集大量沿海滩涂开发工程，河口上游泥沙的淤积，使河流的悬浮物浓度增大；而新安江(千岛湖)段是国家级重点风景名胜区，生态保护较好，土地利用变化较少，使得其悬浮物浓度低(8 mg/L以内)且变化较小。

由于采样日期(2016年12月7日)处于钱塘江低潮位，建立的Chla和TSM模型用于卫星过境时间相近的2016年11月28日的GF-1/WFV数据效果较好。由于钱塘江的潮汐浪涌水流冲击力非常强，给高潮位现场采样带来很大难度，因此现阶段我们的研究主要是围绕钱塘江低潮位情况下展开的，而对于高潮位遥感影像进行测量、建模反演有待今后进一步研究。

5 结 论

(1)基于拟合的GF-1/WFV各波段光谱构建钱塘江流域杭州段叶绿素a和总悬浮物浓度遥感反演模型，其中基于B1-B3差值组合的二次多项式模型可以用于水体叶绿素a浓度的遥感估算，基于B3波段的二次多项式模型可以用于水体总悬浮物浓度的遥感估算，表明WFV数据的波段设置能够较好地用于河流水体水质参数的遥感估算，为高分一号卫星数据遥感监测河流水色参数做出了有益探索，也为后续基于高分遥感影像的河流水质监测的理论和方法研究提供了重要的参考依据。

(2)基于GF-1/WFV数据进行钱塘江流域杭州段叶绿素a和总悬浮物浓度分布制图能够较好地显示其空间分布差异，河口区域总悬浮物浓度和叶绿素a浓度较高，内河段或千岛湖区域水质较好，GF-1/WFV的高空间分辨率对钱塘江河流水体的水质变化细节特征描述得较为详细。

(3)由于钱塘江恶劣的潮汐浪涌水动力特点，为现场采样测量带来了困难，因此对小数据量样本建模的代表性及准确性进行验证是下一阶段研究的重点。另外通过FLAASH的GF-1/WFV大气校正对叶绿素a和总悬浮物浓度分布变化进行精确定量反演与验证也有待进一步深入研究。

致谢：野外试验观测得到浙江省水文局相关人员的支持和帮助，室内试验分析得到国家海洋局第二海洋研究所卫星海洋环境动力学国家重点实验室相关人员的帮助，在此表示感谢。