基于QAA算法的近岸和内陆水域固有光学量反演

Inherent Optical Properties Inversion Coastal and Inland Waters Based on QAA Algorithm

摘要: 水资源在人类的生活环境中起着十分重要的作用，固有光学特性(IOPs)作为反映水生生态环境的重要因素可以使用QAAs算法进行估算。然而，大多数算法依赖412 nm波段估算水组分吸收系数。新半分析算法(QAA-LN)基于遥感反射率比值估算内陆和近岸水体的IOPs。该算法主要由两个部分组成。首先，选择常见的665 nm波段作为参考波段反演总吸收系数a和后向散射系数b_bp。其次，通过490 nm和665 nm波段遥感反射率比值计算浮游植物吸收系数(aph)和颗粒和有色溶解有机物(adg)。使用收集的近岸和内陆水域数据集对模型进行评估，a (λ)、b_bp (λ)、a_ph (λ)和a_dg (λ)的平均绝对误差(MAE)分别为0.20 m⁻¹、0.18 m⁻¹、0.19 m⁻¹、0.28 m⁻¹，结果表现良好。还将新算法用于实测数据匹配的Landsat-8 OLI卫星数据，也能展现出较好的性能，并应用到澳大利亚Hume水库。QAA-LN有潜力作为一种简单有效的算法，未来将扩展到其他卫星监测全球内陆和近岸水域IOPs的变化。

Abstract: Water resources play a crucial role in human living environments, and inherent optical properties (IOPs) serve as important indicators of aquatic ecological conditions that can be estimated using QAA algorithms. However, most algorithms rely on the absorption coefficients of water constituents estimated at the 412 nm wavelength band. The new semi-analytical algorithm (QAA-LN) estimates IOPs for inland and coastal waters based on the ratio of remote sensing reflectance. This method consists of two main components. Firstly, it selects the commonly used 665 nm band as a reference to retrieve the total absorption coefficient (a) and the backscattering coefficient (b_bp). Secondly, it calculates phytoplankton absorption coefficient (a_ph) and the absorption by particles and colored dissolved organic matter (a_dg) using the ratio of remote sensing reflectance at 490 nm and 665 nm. The model was evaluated using the collected coastal and inland water data sets. The mean absolute error (MAE) of a (λ), b_bp (λ), a_ph (λ) and a_dg (λ) were 0.20 m⁻¹, 0.18 m⁻¹, 0.19 m⁻¹, and 0.28 m⁻¹, indicating satisfactory performance. The new algorithm was also applied to the Landsat 8 OLI satellite data in conjunction with the measured data, demonstrating favorable performance, particularly in its application to the Hume Lake in Australia. QAA-LN has the potential to serve as a simple and effective algorithm, with plans for future expansion to monitor changes in IOPs of global inland and coastal waters using other satellites.

1. 引言

固有光学特性(IOPs)是指海洋和内陆水体、浮游植物、有色溶解物、非藻类颗粒物的吸收和散射特性[1]，与水体组分浓度及其光学特性相关。固有光学量可用于估算叶绿素浓度[2]、悬浮物浓度[3]、浮游植物生物量及初级生产力[4]、碳储量[5]和模拟其他生物地球化学过程[6]，对于海洋和内陆水体监测至关重要。通过野外测量和遥感技术准确获取固有光学量，是水色研究的基础[7]。

目前国际上已有多种水体固有光学量遥感反演算法。其中，准分析算法(Quasi-analytical Algorithm, QAA)是水色遥感中应用最广泛的半分析算法，已在国内多个水域得到了成功应用[8]-[10]。该算法主要由两部分组成，首先反演总吸收系数(total absorption coefficients, a)和后向散射系数(backscatter coefficient, b_bp)，其次将总吸收系数分解为水体各组分的吸收系数，包括浮游植物吸收系数(absorption coefficient of phytoplankton, a_ph)和颗粒和有色溶解有机物吸收系数(absorption coefficients of NAP and colored dissolved organic matter (CDOM), a_dg)。

近年来，QAA算法在参考波段[11]-[14]和经验参数[15]-[17]方面经历了几次改进，其性能得到了较高提升。但仍存在以下两方面的问题：(1) 在光学性质复杂的近岸和内陆水域中，受水域环境参数影响，水体各组分吸收系数反演效果存在较大不确定性[18]-[23]。(2) QAA算法反演过程中需要410 nm波段，而大多数卫星传感器缺少该波段。因此，如何在缺少410 nm波段的前提下在复杂水域环境中准确反演水体组分吸收系数是半分析算法面临的难题。

Landsat-8卫星搭载的陆地成像仪(Operational Land Imager, OLI)和能够提供涵盖可见光至近红外的9个波段，具有较高的空间分辨率(30 m)和光谱分辨率，满足了对沿海和内陆水域的调查且在对地观测研究中得到了广泛应用[24]-[26]。长时间数据积累和高分辨率，使得OLI传感器更适合近岸和内陆水体参数的遥感反演[26] [27]。

因此，本文提出一种适合于实测和遥感数据的半分析算法(QAA-LN)，试图反演近岸和内陆水域固有光学量。本研究的主要目的是：(1) 构建一个新半分析算法用于反演固有光学量。(2) 利用全球近岸、内陆水域数据集和Landsat-8 OLI数据，评估新算法的性能。(3) 基于2016~2023年的Landsat-8 OLI数据，分析澳大利亚Hume水库IOPs时空变化特征。

2. 数据

2.1. 实测数据

Table 1. Dataset name, measurement parameters, and sample size

表1. 数据名称、测量要素和样本数量

数据	要素	样本数量
PACE数据集	Rrs (λ), a (λ), aph (λ), ag (λ),	756
	ap (λ), ad (λ), bb (λ)
澳大利亚水域数据集	Rrs (λ), ag (λ)	117
	aph (λ), ad (λ), bb (λ),
渤海数据集	Rrs (λ), a (λ), Zsd (λ), Kd (λ)	75

Table 2. Definitions of symbols used in this article

表2. 本文所用的符号定义

名称	定义
Rrs (λ)	水上遥感反射率
rrs (λ)	水下遥感反射率
a (λ)	总吸收系数，aw (λ) + anw (λ)
aw (λ)	纯水吸收系数
anw (λ)	非水体组分吸收系数，ad + aph + ag
ad (λ)	非藻类颗粒物吸收系数
aph (λ)	浮游植物吸收系数
ag (λ)	有色溶解有机物吸收系数
adg (λ)	颗粒和有色溶解有机物吸收系数，ag (λ)+ ad (λ)
ap (λ)	颗粒物吸收系数，ad + aph
bbw (λ)	纯水后向散射系数
bb (λ)	总后向散射系数，bbw (λ)+ bbp (λ)
bbp (λ)	悬浮颗粒后向散射系数
Ybbp (λ)	悬浮颗粒后向散射系数光谱斜率
S	有色溶解有机物光谱斜率
Chl-a	叶绿素a
Zsd	透明度
Kd	漫射衰减系数

本文利用Casey等[28]提供的PACE实测数据、Qing [29]等提供的渤海海域实测数据和Drayson等[30]提供的澳大利亚水域实测数据等3个独立数据集，评估了QAA-LN算法。其中PACE数据集主要分布在全球近岸海域，渤海海域数据集主要分布在中国渤海近岸海域，澳大利亚水域数据集主要对内陆湖泊、河口、河流等测量。每组数据均由实测R_rs和IOPs数据组成，各参数定义和样本数量见表1和表2。

PACE数据集收集了2004年10~12月和2009~2016年旧金山湾、苏必利尔湖、佛罗里达河口及太平洋海域756组实测数据(见表1)。数据获取及数据详细说明见[28]。澳大利亚水域数据集提供了2013~2021年澳大利亚34个湖泊、河流、大坝、河口等内陆水体的316组生物光学参数[30]，经质量控制最终选择117组数据(见表1)用于模型验证。渤海数据集由Qing [29]等提供，该数据集收集了2005年6月和9月渤海海域75组数据(见表1)，所测参数包括R_rs、Z_sd、a和K_d。三组数据集的数据处理方法近似，除部分可直接获取实测值外，均采用紫外分光光度计测量总颗粒物、有色溶解有机物和非藻类颗粒物的吸收度。对总颗粒物样品甲醇萃取去除浮游植物的影响，将萃取后的空白样本和样本分别置于积分球上、测量获得非藻类颗粒物光谱吸收系数。将冷藏样本升温至室温并置于10% HCL中浸泡、测量得到有色溶解有机物的吸收系数。

2.2. Landsat-8 OLI数据

本文所用卫星数据为2016~2023年澳大利亚4个内陆湖泊、河流、水库区域的Landsat-8 OLI Level-1C数据，该数据由USGS提供(https://earthexplorer.usgs.gov)。经质量筛选，共获取114幅无云影像(见表3)。利用ACOLITE模型对Landsat-8 OLI数据进行大气校正。该模型默认使用暗光谱拟合算法对暗像元处理，减少传感器、气溶胶等其他因素等误差来提高大气校正精度[31]。利用与卫星数据匹配的19组实测R_rs评估大气校正结果，发现大气校正后的R_rs与实测R_rs之间存在较好的一致性(R² = 0.22，MAE = 0.02 m⁻¹，RMSE = 0.02 m⁻¹，见图1)，表明了ACOLITE模型的可靠性。

Figure 1. The average R_rs of Landsat-8 OLI is compared with the measured data in Australia

图1. 澳大利亚实测数据与Landsat-8 OLI的平均R_rs比较

Table 3. Landsat-8 OLI imagery of Australian waters matching measured data

表3. 澳大利亚水域匹配实测数据的Landsat-8 OLI影像

日期	数据
Landsat-8 OLI	2017-03-07	2017-04-28	2017-11-14	2020-02-23
	2020-07-17

3. 方法

3.1. QAA-LN模型

QAA-LN算法流程可分为两个部分(见图2)。第一部分参考QAA-V6算法推导a (λ)和b_bp (λ)，设置参考波段为665 nm；第二部分使用经验公式，通过a_g计算a_dg和a_ph。

Figure 2. QAA-LN algorithmic framework, in which red, yellow, blue, and purple represent input parameters, empirical formulas, analysis formulas, and output parameters, respectively

图2. QAA-LN算法流程图，其中红色、黄色、蓝色和紫色分别表示输入参数、经验公式、分析公式和输出参数

3.2. 精度评价

使用均方根误差(RMSE)、平均绝对误差(MAE)和决定系数(R²)评估算法精度，计算公式(1) (2) (3)如下：

$\text{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(y_i-x_i\right)}^2}}{n}}$ (1)

$\text{MAE}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n\left|y_i-x_i\right|}}{n}}$ (2)

$R^2=\frac{{\left[{\displaystyle {\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}\right]}^2}{{\displaystyle {\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2{\displaystyle {\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}$ (3)

其中x_i为实测值，y_i为反演值，n为样本数。

4. 结果

4.1. 最佳敏感波段确定

参考Kowalczuk [32]，基于实测数据，选择443 nm、490 nm、560 nm、665 nm波段，分析了a_g与不同波段比值之间的敏感性(见图3)，构建了估算a_g的经验公式。结果表明，a_g与R_rs (665)/R_rs (490)之间的相关性最高，R²约0.76。a_g估算公式如下：

$a_g={10}^{c_1*\log \left(R_{rs}\left(665\right)/R_{rs}\left(490\right)+c_2\right)}$ (4)

其中C₁ = 0.6978，C₂ = 0.0077。

Figure 3. Sensitivity of the ratio of bands

图3. 波段比值敏感性

4.2. 基于实测数据的算法精度评估

基于PACE、渤海、澳大利亚等三个实测数据集，对改进的QAA-LN算法进行了精度评价。为了将该算法应用于Landsat-8 OLI数据，评估了443 nm、490 nm、560 nm和665 nm波段处的反演精度。

结果表明，QAA-LN算法在上述四个波段处反演精度高(见图4、见图5)。对于a而言，在560 nm和665 nm处整体反演精度较高(R² = 0.11~0.70、MAE = 0.10~0.11 m⁻¹、RMSE = 0.21~0.31 m⁻¹)。b_bp基本分布在1:1线两侧，其中R²、MAE和RMSE精度值最高分别为0.90、0.04 m⁻¹和0.09 m⁻¹。a_dg的反演结果与实测值相比表现为略微高估，不同波段处R²和RMSE精度值近似分别在0.56~0.65和0.08~1.27 m⁻¹间，MAE精度最高值为0.05 m⁻¹。与a_dg反演结果相反，a_ph出现明显低估现象，尤其是443~560 nm波段。与上述三个波段相比精度有所降低，但MAE精度值也在0.14 m⁻¹~0.26 m⁻¹之间。上述结果表明了该算法的可靠性。

Figure 4. Comparison of measured values and inversion values of a and b_bp

图4. a和b_bp实测值与反演值比较

Figure 5. Comparison of measured values and inversion values of a_dg and a_ph

图5. a_dg和a_ph实测值与反演值比较

4.3. 卫星反演精度评估

将上述算法应用于大气校正后的Landsat-8 OLI数据，得到基于卫星数据的反演结果。将反演结果与准同步的实测数据进行对比，发现a的反演效果更好，数据大多数分布在1:1线上，与实测数据基本保持一致(见图6、见图7)。其中，a (665)反演精度最高，MAE和RMSE值均低于0.32 m⁻¹，R²达到0.25。b_bp整体表现为略微高估，R²最低也达到0.57且MAE和RMSE精度值均在0.03 m⁻¹左右。a_dg和a_ph反演结果分别表现为低估和高估，算法均在665 nm处反演精度最好，其MAE、RMSE精度值最高可达到0.08 m⁻¹~0.12 m⁻¹和0.09 m⁻¹~ 0.17 m⁻¹。因此，将算法应用于遥感数据，具有一定的可靠性。

Figure 6. Comparison a and b_bp remote sensing inversion values and measured values

图6. a和b_bp遥感反演值与实测值比较

Figure 7. Comparison a_ph and a_dg remote sensing inversion values and measured values

图7. a_ph和a_dg遥感反演值与实测值比较

4.4. 遥感应用

以澳大利亚Hume水库为例，将算法应用于2016~2023年Landsat-8 OLI遥感影像，分析其年均a_ph (665)时空分布特征。总体而言，Hume水库IOPs表现出显著的时空差异(见图8、见图9)。从时间上看，整体主要分为2个阶段。其中2016年~2017年、2018年~2019年和2021年~2023年a_ph (665)值呈下降趋势，可能是人类活动减少导致。2019年~2021年呈上升趋势，在2021年出现最大值。参考吴月圆等[33]和孙军等[34]结果，可能是由于当年发生火灾，气温升高促进了浮游藻类的生长。从空间上看，a_ph (665)最大值出现在东北部分可能是由于离岸较近，人类活动较多导致的IOPs变化。异常值也常出现在东北部水域。a_ph (665)最小值分布在水库南部。此外，我们还发现在水库沿岸地区变化明显原因是浅水处湖底的反射率影响了算法精度[35] [36]。

Figure 8. Time variation of IOPs in Hume Lake from 2016 to 2023

图8. 2016~2023年Hume水库IOPs时间变化

Figure 9. Spatial variation of IOPs in Hume Lake from 2016 to 2023

图9. 2016~2023年Hume水库IOPs空间变化

5. 讨论

5.1. 不同算法对比

Figure 10. Comparing the measured data of a with the inversion data, the inversion results of QAA-LN, GIOP and GSM algorithms are from left to right

图10. a的实测数据与反演数据比较，从左到右依次为QAA-LN、GIOP和GSM算法反演结果

Figure 11. Comparing the measured data of bbp with the inversion data, the inversion results of QAA-LN, GIOP and GSM algorithms are from left to right b_bp

图11. b_bp的实测数据与反演数据比较，从左到右依次为QAA-LN、GIOP和GSM算法反演结果

Figure 12. Comparing the measured data of a_ph with the inversion data, the inversion results of QAA-LN, GIOP and GSM algorithms are from left to right

图12. a_ph的实测数据与反演数据比较，从左到右依次为QAA-LN、GIOP和GSM算法反演结果

Figure 13. Comparing the measured data of a_dg with the inversion data, the inversion results of QAA-LN, GIOP and GSM algorithms are from left to right

图13. a_dg的实测数据与反演数据比较，从左到右依次为QAA-LN、GIOP和GSM算法反演结果

在近岸和内陆水域进行，对比了QAA-LN、GIOP [37]和GSM [38]算法性能(见图10、见图11、见图12、见图13)，发现QAA-LN性能最优，其次为GIOP和GSM。具体为，在560 nm和665 nm波段处，QAA-LN算法反演的a和b_bp与实测数据基本保持一致(R²、MAE和RMSE取值范围分别为0.18~0.90、0.04~0.10 m⁻¹和0.09~0.21 m⁻¹)，a_dg和a_ph反演结果分别表现为低估(R²、MAE和RMSE取值在0.05、0.15 m⁻¹、0.27 m⁻¹左右)和高估(R²、MAE和RMSE取值范围分别为0.56~0.66、0.05~0.19 m⁻¹和0.08~0.30 m⁻¹)。GIOP算法反演精度低于QAA-LN，R²、MAE和RMSE取值范围分为别0.01~0.74、0.08~1.27 m⁻¹和0.07~0.77 m⁻¹。GSM算法精度最低，存在高估现象。

5.2. 误差传播

本研究讨论了可能导致IOPs反演精度不确定的因素。总的来说，参考波段的选择(λ0)、参考波段处吸收系数(a (λ0))经验公式的建立、Y_bbp以及S都会影响吸收系数的反演精度。具体而言，对于a来说，a (λ0)的计算精度至关重要，需要根据实测数据建立经验模型[39]。在不同水体类型中，不同参考波段会影响bb (λ)精度，尤其是当参考波长转移到较长波段时(708 nm、710 nm)，总吸收系数以纯水吸收系数为主，非水分子吸收系数的影响可忽略不计，因此无需经验公式(Δa (λ0))即可获得总吸收系数(a (λ0) ≈a_w (λ0)) [29] [40] [41]。Y_bbp在推导IOPs具有二阶重要性[39]。在计算太湖水域[42] [43]的Y_bbp值时发现实测值与推导值相差较大。Y_bbp对b_bp (λ)的影响受(λ0/λ) Y_bbp的组合变化。Yang [44]在其研究中对其进行了详细的测试，发现当Y_bbp实测值大于推导值，反演的b_bp与反演值相比偏低；当Y_bbp实测值小于推导值，b_bp将比实测高。在此处产生的错误将会进一步影响a的精度。因此，如要想将算法应用到近岸和内陆水域，还需要对Y_bbp值进行准确估算。此外，S也会随着公式将误差传递到a_dg和a_ph。在以往的研究中[37]，S通常使用固定值，他们认为模型中函数截距可看作S值。后续研究[45]表明S在0.01~0.02 nm间变化，因此在后来的QAA模型中，使用0.015作为平均值并考虑了红绿波段比提高精度。不同水域的光谱斜率存在空间差异，当S发生变化时，会引起结果的额外误差[46]-[49]。因此，如何在复杂环境中保证算法具有良好精度也是当前算法改进需要注意的难题。

6. 总结

本文基于全球近岸、内陆水域数据集，提出了一种基于QAA算法的近岸和内陆水域半分析模型(QAA-LN)。该算法以中高分辨率卫星通用波段665 nm作为参考波段，基于实测遥感反射率R_rs (665)和R_rs (490)，计算了a_dg和a_ph，提高了近岸和内陆水域IOPs反演精度。基于实测数据的独立验证结果表明，该模型实测数据反演精度：R² = 0.60、MAE = 0.24 m⁻¹、RMSE = 0.63 m⁻¹。卫星数据反演精度：R² = 0.06、MAE = 0.60 m⁻¹、RMSE = 1.57 m⁻¹。与GIOP、GSM半分析算法相比，该算法在估算总吸收系数和浮游植物吸收系数表现出良好性能。将算法应用于Landsat-8 OLI卫星数据，能有效地反演澳大利亚Hume水库IOPs的空间分布特征。

综上所述，该方法能够较好地反演近岸和内陆水域IOPs。然而，由于光学特性的多变性，在优化模型时需要更多的数据和工作寻找变量限制进一步优化算法。模型也可应用于与Landsat-8 OLI相似的高分辨率卫星传感器Sentinel-2MSI等。本研究为近岸和内陆水域监测和管理提供了一种有用的方法，结合Sentinel-2 MSI卫星数据得出的IOPs产品将成为探索复杂光学特性的近岸和内陆水域水组分长时间序列变化的趋势。

致 谢

感谢导师青松和阿如娜博士为本研究提供了宝贵的研究思路和论文撰写过程中给予耐心的指导。感谢Casey K A、Drayson N等提供的实测数据。感谢美国地质调查局(USGS)提供Landsat-8 OLI数据。

基金项目

本研究获得国家自然科学基金(41961057)、内蒙古自治区应用数据中心自主研究重点项目(ZZYJZD2022001)和内蒙古自治区高等学校青年科技英才支持计划项目(NJYT23016)项目支持。