资源一号02D卫星高光谱数据水体透明度反演研究

周毅 刘瑶 田淑芳

(1 中国地质大学(北京) 地球科学与资源学院，北京 100083)(2 自然资源部国土卫星遥感应用中心，北京 100048)

透明度作为反映水体清澈/浑浊程度的水质参数，被广泛应用于湖泊学和海洋学研究中，对评估水体内水生植物多样性、生产力和水体营养程度等方面具有重要意义[1]。过去几十年来，受人口增长、经济发展中不合理的开发利用及污染等因素影响，我国许多内陆水体(湖泊、河流和水库)呈富营养化状态，水体质量急剧下降，带来水体生态系统的结构破坏与功能退化，对人民生活及国家经济造成巨大损失和不良影响[2]。近年来，国家为改善水生态环境，开展了重点流域保护修复、重点工程水质保障等一系列工作。有必要定期开展内陆水体透明度等水质参数的监测工作，服务于生态修复措施的水环境改善效果评价，以及为后续水资源的保护和管理提供决策数据支持。透明度通常定义为塞氏盘(Secchi Disk)在水中垂直下放至刚刚在视野中消失的深度。由于塞氏盘获取成本低且操作简便，使用塞氏盘进行水体透明度测量已有上百年的历史。然而，通过布设地面监测点或者船载测量的方式难以实现对大区域和多湖泊的透明度重复观测。相比而言，能够周期性稳定获取和实现空间覆盖的遥感卫星，成为了开展长期大规模水体透明度反演的重要数据源。

国内外已有众多学者基于遥感数据开展水体透明度反演研究，常见方法有经验/半经验模型和半分析模型。文献[3]基于陆地卫星-5(Landsat-5)专题成像仪(TM)数据的蓝波段、红蓝波段反射率比值，使用相关系数矩阵、多元逐步回归等方法完成了1973～1998年间美国大型湖泊的透明度评价研究；文献[4]使用Landsat-5和陆地卫星-7(Landsat-7)增强型专题成像仪(ETM+)数据的蓝红波段，通过建立经验模型对1990～2010年间美国缅因州湖泊开展透明度反演研究；文献[5-6]建立了一种基于半分析算法(Quasi Analytical Algorithm，QAA)的水体透明度反演模型，可应用在包括中分辨率成像光谱仪(MODIS)、中等分辨率成像频谱仪(MERIS)、陆地卫星-8(Landsat-8)等卫星的多种传感器上实现水体透明度的反演。国内方面，文献[7]使用Landsat-5中红蓝波段建立自然对数与透明度自然对数的线性模型，对2004年6个时期鄱阳湖的水体开展透明度反演，取得了较好精度；文献[8]应用Landsat-8 陆地成像仪(OLI)数据的红、绿波段建立对数模型，反演了2013年的东平湖水体透明度，取得了19.77%的总体相对误差；文献[9]基于透明度和水体漫衰减系数关系建立半分析模型，并通过MODIS、近海高光谱成像仪(HICO)数据对近岸水体透明度进行反演。

目前，水体透明度反演研究多以多光谱数据为数据源，采用星载高光谱数据开展研究相对较少。相较于多光谱数据，高光谱数据具有更高的光谱分辨率，在用于水体透明度反演时，能够提供更丰富的波段选择，且具备更精准的光谱特征，因此在水质参数反演精度的提升方面具有优势和潜力。

我国于2018年5月9日成功发射高分5号(GF-5)卫星，其搭载新一代高光谱相机(Advanced Hyper Spectral Imager，AHSI)的光谱覆盖范围为0.40～2.50 μm，空间分辨率为30 m，可见光近红外(VNIR)与短波红外(SWIR)分别有150和180个波段，光谱分辨率分别为5 nm和10 nm。此外，资源一号02D卫星(又称为5米光学业务卫星)也于2019年9月12日成功升空。该星搭载的AHSI，具有与GF-5卫星相同的光谱覆盖范围和空间分辨率，可见光近红外波段有76个通道(光谱分辨率为10 nm)，短波红外波段包括90个通道(光谱分辨率为20 nm)。我国后续还将发射搭载高光谱相机的卫星，通过多星组网，未来有望能够基于星载高光谱数据开展内陆水体水质监测。

本文以北京市官厅水库为研究区，以资源一号02D卫星高光谱影像为数据源，应用半经验和半分析模型方法反演研究区的水体透明度。通过对结果精度进行评价，分析资源一号02D卫星高光谱数据在反演水体透明度方面的表现与应用潜力。

1 材料与方法

1.1 研究区概况

官厅水库位于北京西北部永定河上游(见图1)，水库面积约130 km2。1970年以来，受入库流量下降、污染负荷增加等影响，官厅水库水质不断恶化，于1997年退出北京饮用水水源地之列；后经治理，水库在2007年再度成为北京备用饮用水水源地。近年来，随永定河综合治理和生态修复推进，官厅水库水质逐年好转[10]。

图1 研究区示意图Fig.1 Map of study area

1.2 数据及处理

2020年8月13日在官厅水库开展了与资源一号02D卫星同步的地表试验，获取了地表实地的光谱数据和水质参数(叶绿素a、悬浮物浓度、透明度等)。高光谱遥感数据的成像时间为2020年8月13日上午11时21分，地面试验则在当天10～14时内完成，以保证星地数据获取的一致性。本次试验在研究区范围内布设了18个均匀分布的点位(见图2)，其中实测水体透明度范围为0.4～0.7 m，均值为0.58 m。

图2 水体透明度实测位置图Fig.2 Location map of measured water transparency

对获取的星载高光谱数据，进行了辐射定标、大气校正、水体范围提取和遥感反射率(Rrs)计算等一系列预处理工作，具体的处理过程如图3所示。

图3 高光谱数据处理过程Fig.3 Hyperspectral data process flow

1)辐射定标

通过ENVI5.3软件，将原始影像的VNIR数据(76个波段)和SWIR数据(90个波段)合并为一个文件，应用Apply Gain and Offset工具导入绝对定标参数后进行辐射定标。

2)大气校正

FLAASH大气校正模型是目前在多/高光谱影像中使用效果最好的大气校正工具之一，对于地物波谱信息可以高保真地恢复，简洁快速地获取比较准确的物理模型参数[11]。因此本文在ENVI5.3中应用FLAASH模块对影像数据进行大气校正。综合考虑数据特点、成像时间及研究区位置，本文在FLAASH模块中选择中纬度夏季(Mid-Latitude Summer)大气模型，水体反演波段设置为940 nm，并选用乡村(Rural)气溶胶模型，以及2-波段(考夫曼-泰伦)(2-Band(K-T))气溶胶反演方法。

3)水体范围提取

本文选择归一化水体指数(Normalized Difference Water Index，NDWI)方法进行水体范围提取[12]，其公式为

(1)

式中：ρGreen、ρNIR分别代表绿波段和近红外波段的地表反射率。

式(1)最初以Landsat TM数据进行计算，而资源一号02D卫星高光谱数据信噪比相对较低，选择单一波段计算容易受到噪声的影响[13]。因此，本文参照Landsat TM数据的绿光(0.52～0.60 μm)和近红外波段范围(0.76～0.90 μm)，通过波段合并，分别计算资源一号02D卫星数据在Landsat TM这两个波段范围内的地表反射率平均值，用于归一化水体指数的计算。在此基础上，运用直方图双峰法(Histogram Bimodal Method，HBM)分离水体和其他类型地物，并将小于1000像元(约0.9 km2)的微小水体消除，只保留主体水域，从而实现影像中水体范围的提取。

4)遥感反射率计算

遥感反射率计算公式为

(2)

式中：λ为波长，ρ(λ)为经过FLAASH大气校正后的地表反射率，ρsky(λ)为天空光反射率，Lu(λ)为水面上行辐亮度，Ld(λ)为水面下行辐亮度，Lsky(λ)为天空光辐亮度。

即使是反射率相对较高的浑浊水体，在短波红外的反射率也趋于0，因此通常认为图像中短波红外的反射率是来自天空光的影响，将短波红外的反射率用于天空光的去除[14]。但是，由于AHSI的辐射定标系数是基于地表亮目标计算得到，而水体反射率相对较低，造成水体像元在短波红外波段噪声相对较大，不适宜用于天空光的纠正。因此，本文采用忽略天空光的影响，以ρ(λ)/π作为遥感反射率的近似值。尽管采用这样近似计算的结果可能略大于实际的遥感反射率，但由于通常认为天空光的影响与波长无关，因此，近似计算的遥感反射率与真实反射率具有一致的光谱形状。基于此，本文采用波段比值的经验模型，和半分析模型进行透明度反演。这两类方法对于透明度的反演，主要取决于光谱的波形特征，受遥感反射率的量级影响较小[15]。

1.3 水体透明度反演方法

本文分别运用半经验模型和半分析模型方法进行官厅水库水体透明度反演。

1.3.1 半经验模型

根据文献[16]研究成果，半经验模型方法主要通过对地面同步点位的影像光谱和地表测量的透明度进行相关性分析，选择合适的波段或波段组合进行建模。本文通过分析资源一号02D卫星图像和地表采集的水体光谱特征，选择中心波长为653 nm的波段，采用对数回归模型(见图4)得到水体透明度反演公式为

Zsd=-0.438×ln (Rrs(653))-1.335 8

(3)

式中：Rrs(653)为中心波长为653 nm波段处的遥感反射率。

图4 对数回归模型Fig.4 Logarithmic regression model

1.3.2 半分析模型

通过半分析模型反演水体透明度主要步骤为：①估算固有光学量(IOPs)，即总吸收系数a和后向散射系数bb；②基于a和bb反演漫衰减系数(Kd，m-1)；③用Kd和Rrs确定Zsd。

本文分别运用第5版及第6版半分析算法(QAAv5、QAAv6)[17]计算总吸收系数a和后向散射系数bb；并采用更新后的半分析模型[18]对漫衰减系数Kd进行反演

m1×(1-m2×e-m3×a(λ))bb(λ)

(4)

式中：θs为太阳天顶角，bbw(λ)为纯水的后向散射系数，m0-3与γ分别取0.005、4.26、0.52、10.8和0.265的固定值，a(λ)为总吸收系数，bb(λ)为后向散射系数。

对于水体透明度的反演，本文通过公式[19]进行计算

(5)

式中：Kd(tr)表示水体在可见光谱范围(410～665 nm)上的漫衰减系数，Rrs(tr)表示该波长下相应的遥感反射率。

1.4 评价方法

本文通过建立线性回归模型来分析反演结果与实测值之间的相关性，并计算绝对误差(AE)、相对误差(RE)以及均方根误差(RMSE)用于对结果的精度评价为

IAE=|dm-dp|

(6)

(7)

(8)

式中：dm和dp分别表示透明度的实测值和反演值；n为样本数。

2 结果与分析

2.1 不同反演结果与实测值对比

本文分别运用1.3节中的方法采用资源一号02D卫星数据对研究区水体透明度进行计算，并提取了与实测点位匹配的相应反演值。应用线性回归模型对反演结果与实测值进行了相关性分析并得到相应的相关系数值(见图5)。

图5 线性回归分析图Fig.5 Linear regression analysis results

此外，运用1.4节中各误差分析方法计算得到各反演结果的误差分析情况(见表1)。

表1 反演结果误差分析表Table 1 Error analysis of inversion results

根据线性回归分析结果可以看出，运用半经验模型方法反演得到的水体透明度结果与实测结果相关性最高，R2为0.473 3；而基于QAAv5与QAAv6的半分析模型方法反演结果相关性较低，其中基于QAAv5计算的反演结果与实测值相关性最低，R2为0.419 3。图5(b)和(c)中，由于QAAv5与QAAv6算法主要是针对一类大洋水体提出，是以550～560 nm作为标准波长；将其应用于浑浊水体，会造成对总吸收系数和后向散射系数的低估，进而低估漫衰减系数，导致反演的水体透明度偏高[20-21]，因此拟合线低于1∶1线。通过表1可知，半经验模型取得的反演结果，相较2种半分析模型方法结果，整体精度较高，绝对误差平均值为0.05 m，平均相对误差为9.5%，而RMSE仅为0.07 m；两种半分析模型得到的反演结果精度相对较低，其中，基于QAAv5计算的反演结果精度较差，最大相对误差达到了174.0%，其平均相对误差为94.7%，RMSE为0.55 m。而基于QAAv6反演结果平均相对误差及RMSE分别为31.4%、0.19 m，相较QAAv5结果精度稍高，主要原因是QAAv6算法对浑浊水体的情况进行了适当的优化，对浑浊的内陆水体的适用性稍有提升。

上述3种方法反演的透明度结果与实测结果的R2都比较低，这主要是由于本研究区的同步点位数量有限；另一方面，反演结果的精度分析可以看到，表明该方法在应用于官厅水库透明度反演中具有较好的效果。

2.2 官厅水库水体透明度分析

通过上述分析，运用半经验模型方法反演得到的水体透明度具有较高精度，因此本文基于该方法得到官厅水库水体透明度(Zsd)状态图(见图6)。

图6 官厅水库Zsd状态图Fig.6 Zsd state of Guanting reservoir

从图6可以看出，官厅水库水体透明度范围在0～1 m之间，相比东北部水体，西南部水体透明度更高，说明该区域水体质量更好。该结果与实地测量得到的水质状态基本一致。

3 结论

本文以北京市官厅水库为研究区，采用资源一号02D卫星高光谱数据，应用半经验模型、半分析模型方法进行研究区的水体透明度反演。通过反演结果与实地测量结果的相关性分析和精度评价，得到主要结论如下：

(1)本文分别使用半经验模型、基于QAAv5与QAAv6的半分析模型反演得到水体透明度，通过结果与实测数据相关性及精度评价，半经验模型方法得到的反演值与实测值相关性最高，相关系数为0.473 3，且结果精度也最高，平均绝对误差、平均相对误差及RMSE分别为0.05 m、9.5%、0.07 m；其次为基于QAAv6半分析模型方法的反演结果，平均绝对误差、平均相对误差及RMSE分别为0.17 m、31.4%、0.19 m；而基于QAAv5的半分析模型反演结果精度较差。

(2)基于资源一号02D卫星的高光谱数据，可实现水体范围的提取并用于内陆浑浊水体的透明度反演。相比较常规的水色卫星，资源一号02D卫星具有更高的空间分辨率，可应用于中小型湖库的水质监测。下一步的研究是验证该数据应用于更多不同光学特性水体中的透明度反演能力。