基于辐射传输模型的巢湖叶绿素a浓度反演

刘文雅， 邓孺孺,2,3， 梁业恒， 吴仪， 刘永明

(1.中山大学地理科学与规划学院，广州 510275； 2.广东省水环境遥感监测工程技术研究中心，广州 510275； 3.广东省城市化与地理环境空间模拟重点实验室，广州 510275)

0 引言

近年来，工农业等污水过度排放导致藻类等水生植物频频暴发，国内乃至全球湖泊富营养化问题亟待解决。湖泊富营养化主要表现为水华暴发，浮游植物呼吸及分解消耗了大量水中的溶解氧，导致水体生物死亡率上升，生态系统失衡。此外，其主要种群微囊藻对人体具有致癌作用[1]，毒素通过饮用水或湖内鱼类威胁人体健康。巢湖为我国五大淡水湖之一，地处于安徽省合肥市东南部，流域内共有33条河流，贯穿7大水系。巢湖流域不仅承载着江淮地区的水循环重担，还维系着当地文化及旅游产业的发展。

近年来巢湖富营养化严重，利用遥感方法的巢湖水质研究始终处于热点。李素菊等[2]于2002年利用高光谱数据，建立反射率比值模型以及一阶微分模型应用于巢湖水体叶绿素a浓度估算，相关系数分别为0.70和0.66； 荀尚培等[3]于2009年利用MODIS最佳通道组合，定义叶绿素a指数IChla反演了巢湖叶绿素a浓度； 杨煜等[4]在2010年通过星地同步实验，基于环境一号卫星高光谱数据建立了3波段反演模型，决定系数为0.868 8； 谢杰等[5]于2010年利用TM数据结合波段组合法，反演并分析6 a间巢湖叶绿素a浓度演变趋势； 陈静等[6]于2012年基于环境一号卫星CCD数据建立模型反演并分析了巢湖叶绿素a季节分布特征； 殷守敬等[7]于2018年将地面观测数据空间插值与遥感反演结合，以高精度的叶绿素a浓度为指标评价了巢湖富营养化空间分布情况。以上研究大多基于实测光谱和浓度的统计关系，虽简单易行，但无法同时反演多个时相的叶绿素a浓度，普适性较低。辐射传输模型是基于水质因子的光学特征，直接建立行星反射率与物质浓度的物理方程，其演算过程完全依托影像，摆脱了实测数据的限制。但因模型参数较多，演算复杂，鲜有国内学者应用。而国外已有较多相关成果： Gilerson等[8]于2010年利用水质因子的固有光学性质结合反射光谱对美国内布拉斯加州的几个湖泊进行了叶绿素a浓度提取； Matthews[9]于2011年分析了叶绿素a等若干参数的固有光学特性，确定可用于内陆水域遥感模型的波段； Mishra等[10]于2013年建立了确定光学活性成分吸收系数的分析算法，估算美国密西西比州的若干养殖池塘中的浮游植物浓度； Stumpf等[11]于2016年总结了国外解析算法、半分析算法以及衍生算法估算藻华浓度的优劣，指出辐射传输模型具有更强的稳健性和普适性。故本文依据辐射传输机理，建立叶绿素a浓度反演物理模型； 采用巢湖不同时相的遥感影像验证模型的普适性； 进一步应用模型分析巢湖2016年藻类的时空分布特征。

1 叶绿素遥感反演模型

1.1 叶绿素参与的辐射传输过程

水域上空传感器接收水体辐射的过程如图1所示，到达水面的电磁波一部分发生镜面反射，一部分进入水体，与水分子及各类水质因子发生吸收散射作用，形成不同波段水体光谱的差异[12]。到达水底的部分经过反射回到水面。水面镜面反射、水体散射和水底反射部分的电磁波共同组成水体出水反射辐射[13]。

图1 电磁波与水体及大气相互作用

故传感器接收到的出水辐亮度L(λ)为

L(λ)=Lw+Ls+Lp，

(1)

式中：Lw为整层水体离水辐亮度；Ls为水面镜面辐亮度；Lp为大气的复合辐亮度。

遥感影像预处理可以消除Lp及Ls，处理后水体反射率仅与Lw相关。Lw又包含水体散射部分Lws以及水底反射部分[14]Lg，即

Lw=Lws+Lg。

(2)

水体光谱特征决定于水中组分的吸收散射行为特征，分别以吸收和散射系数表征。对巢湖而言，水体反射率主要影响因子为叶绿素和悬浮泥沙2种。同一物质在同一波段的吸收散射系数为定值，因此可以建立上述2种物质的反射率与光学参数、物质浓度的关系模型，进而选取敏感波段联立方程求解叶绿素a浓度。

1.2 叶绿素a浓度遥感反演模型

根据1.1节中论述，需分析电磁波透过水面后在整层水体中的出水散射以及水底反射2部分。电磁波透过水面进行下行辐射，在水深为h处，厚度为dh的薄层水散射作用可表示为

(3)

(4)

β=βw+Dsβs+Duβu，

(5)

α=αw+Dsαs+Duαu，

(6)

式中：p(Θ)为水层垂直向上的散射比;Θ为水层垂直向上的散射角；E为水深为h处的下行辐照度；E0为太阳入射水面的初始辐照度；θ′为光穿透水面后的入射角；α和β分别为各因子总吸收和散射系数；Ds和Du分别为悬浮泥沙浓度和叶绿素a浓度；αw，αs和αu则分别对应纯水、悬浮泥沙以及叶绿素a吸收系数；βw，βs和βu则分别对应上述3种分子的散射系数。经过下行散射衰减的入射光在离水之前又经过上层水体进行再次衰减[15]，表达式为

(7)

式中φ为水底反射角。对式(7)积分，并将式(4)代入可求得水体散射部分表达式为

(8)

(9)

式中k为消光系数，即吸收系数与散射系数之和。

设水底为均一的朗伯体，Rg为水底反射率，则入射光在水底的离水辐亮度Lg为

(10)

将式(8)和式(10)代入式(2)，已知反射率R与辐亮度L的转化关系为

(11)

最终得到整层水体反射率表达式为

(12)

各波段对应的光学参数均可通过测量或计算获取，故未知数仅有悬浮泥沙浓度Ds、叶绿素a浓度Du、水深h以及水底反射率Rg。分析光学特性，选择敏感波段组合建立并求解方程组，即得到目标水体叶绿素a浓度。

2 巢湖流域叶绿素a浓度提取

2.1 遥感影像预处理

Landsat8增加了用于监测近岸水体和气溶胶的深蓝波段[16]，空间分辨率为30 m。本文所使用的巢湖多时相遥感影像获取于美国地质调查局官网[17]，影像行列号分别为121/038。获取影像共4景，分别成像于2016年6月23日、7月25日、8月10日以及11月14日。采用PCI GEOMATICA软件对其进行辐射定标、大气校正和水陆分离等一系列预处理。

2.1.1 辐射定标

Landsat8遥感图像metadata的头文件可直接用于建立反射率与DN值间的线性关系。第8—11波段在此不予考虑。2016年7月25日影像对应的定标方程为

R(λ)=0.000 021 9DN-0.109 713 5 。

(13)

2.1.2 大气校正

大气对遥感辐射传输过程的影响主要为吸收和散射作用，包括入射光的消光、路径散射、天空漫射光以及因邻区多次散射面混入像元反射而造成的影像模糊等[18]。采用暗像元法进行大气校正[19]。水域之上传感器接收到的短波辐射由水面辐射和大气辐射组成[20]，即

L(λ)=Lr(λ)Tu(λ)+Lu(λ)，

(14)

式中：Lr(λ)Tu(λ)为水面辐射部分；Lr(λ)为初始水面反射光辐亮度，包含天空散射光和水体散射光2部分；Tu(λ)为上行反射光大气透过率；Lu(λ)为大气散射部分。

Lr(λ)表达式为

(15)

式中：Td(λ)为入射光下行透过率；Ed(λ)为天空光辐照度，因测量时间为夏季，天空光部分仅考虑瑞利散射；θ为入射光天顶角；Rw为水体反射率。

Lu(λ)表达式为

(16)

式中：ω为大气总散射系数；P(Θ)为瑞利散射的散射相函数，垂直观测下Θ=θ。则P(θ)公式为

(17)

设ρu为上行路径散射率，则结合定义及一系列变换，得到真实水体反射率Rw与像元反射率R(λ)间的关系为

(18)

则大气校正计算公式为

Rw=aR(λ)+b，

(19)

(20)

(21)

式中Rvd为暗像元反射率，根据文献[19]，认为ρu与Rvd近似相等，其表达式为

(22)

结合各参数定义依次计算出ω，Tu和Td，进而求出a和b，并基于式(19)进行波段运算，得到大气校正后的影像如图2所示。大气校正清除了因天空光散射而产生的蓝色朦胧感，且各地物色彩与实际相符，大气校正效果显著。

图2 大气校正后影像

2.1.3 水陆分离

基于突出水体信息、排除异物同谱现象的原则，选取波段构建分离指数以进行水陆分离。以Landsat8 OLI为例，水体在B6和B7的反射率几乎为0，水体信息最丰富。水体与植被在蓝光—红光波段的反射率相似，直至近红外波段出现植被反射峰，故B5用于分离水体与植被； 水体与土壤仅在绿光波段的反射率相近，将B3用于排除土壤像元； 建筑物反射率在B6骤然增强，将B6作为水体与建筑物分离的特征波段。综上所述，将B3，B5和B6作为提取水体的特征波段，此外，蓝光波段包含的水体信息较为丰富，故将B2列入以增强水体特征。利用4个特征波段组合运算得到的水体提取指数NDI为

(23)

式中B2，B3，B5和B6分别为第2，3，5和6波段反射率。然后，基于NDI值，选取阈值进行二值化，得到水陆分离结果，如图3所示。图3(a)结果可见，巢湖及周边水域均被提取； 图3(b)局部提取结果对比可知，细小河流也没有被漏分，分离精度较高。

(a) 水陆分离结果 (b) 水陆分离局部对比

2.2 光学参数获取

按照参考文献[21]方法获取αw和βw； 参考文献[22]方法获取βs，αu和βu。由于研究区差异，文献[22]忽略了悬浮泥沙的影响，但巢湖属于典型的浑浊二类水体，悬浮泥沙的强散射严重削弱叶绿素a光谱信息，剔除悬浮泥沙光谱干扰是准确反演叶绿素a的关键[23]。因此，在实际计算中，利用遥感影像特殊像元法对特定波段的悬浮泥沙吸收系数进行计算： 在经过预处理后的红光、绿光和蓝光3波段合成的遥感影像上，选择一浑浊严重、且无绿色及黑色物质的水域，代表悬浮泥沙浑浊严重，无有机污染物且叶绿素a浓度较低的特征像元。读取该像元对应Landsat8 OLIB4和B5反射率，分别与已获取的上述5个系数代入式(12)，解得2个波段对应的悬浮泥沙吸收系数αs。

所需光学参数值如表1所示。

表1 水质因子光学参数值

2.3 特征波段选取

利用敏感波段区分叶绿素a及水中其组分是反演的关键[24]。叶绿素a发生光合作用时，在440 nm和670 nm出现吸收峰，750 nm附近出现反射峰[25]，可作为叶绿素a异于其他组分光学特征波段的选取依据。特征波段选取遵循波段间相关性小、排除异物同谱现象的原则[26]。B6和B7这2个波段的水体反射率几乎为0，故排除。叶绿素a在深蓝光、蓝光和红光波段共同表现为强吸收弱散射，故首先将绿光或近红外波段列入敏感波段。绿光波段对应叶绿素a和悬浮泥沙均表现为弱吸收、弱散射，将其排除。考虑到深蓝及蓝光波段受大气散射影响严重，将红光波段列为第2个特征波段。综上所述，最终确定叶绿素a特征波段为B4和B5。

2.4 巢湖叶绿素a浓度反演

巢湖水较深，加之悬浮泥沙影响，穿透水体至水底的入射光几乎为0。近似认为式(12)中h趋近于正无穷大，化简得到最终巢湖叶绿素a浓度计算模型为

(24)

模型未知量包括Du和Ds，基于B4和B5对应的光学参数值，联立并求解方程组，代入预处理后的B4和B5反射率即得整幅影像叶绿素a浓度Du。

3 模型普适性验证

叶绿素a浓度遥感反演算法的普适性强意味着在光学差异较大的水体中皆能取得较高的反演精度[27]。内陆湖不同时相的水体具有较大的光学差异，将该模型同时应用于巢湖不同时相的遥感数据，若能获取较准确的反演结果，可证明模型普适性较强。

3.1 验证数据获取与处理

2组实测浓度数据引用于参考文献[28]和[29]，分别为巢湖2006年7月30日的12个实测点以及2009年3月27日的23个实测点对应的叶绿素a浓度。

同步遥感影像数据来源于地理空间数据云平台(http: //www.gscloud.cn)。受影像含云量限制，2009年3月27日的同步遥感数据使用4月9日的Landsat7 ETM+影像代替。由于受到Landsat7 ETM+机载扫描行校正器故障影响，该影像出现数据条带丢失，使用ENVI软件插件对影像进行去条带处理[30]； 辐射定标利用ENVI软件的Calibration工具完成，其余反演步骤与第2章一致。

3.2 精度及普适性分析

2006年7月30日与2009年3月27日的反演结果如图4所示，2幅影像模型值与实测值对比如图5所示，定量对比结果如表2所示。

(a) 2006年7月30日 (b) 2009年3月27日

(a) 2006年7月30日 (b) 2009年3月27日

表2 模型反演值与实测值数值统计对比

采用统计学方法作为检验模型精度的评判标准[5]，统计标准引用3个指标，即决定系数R2、相对误差RE和均方根误差RMSE，误差结果如表3所示。

表3 模型反演值与实测值的误差

2次模型验证的决定系数R2分别为0.877 8和0.848 8，精度较优，说明巢湖作为光学特性复杂的富营养化浑浊水体，根据统计关系建立的经验及半分析模型强化了除悬浮泥沙外其他组分的影响[5]，增大了估算误差，而辐射传输模型针对内陆湖高浑浊特性，基于光学机理建立叶绿素a与悬浮泥沙浓度间的关系，准确剔除了悬浮泥沙干扰，最终实现叶绿素a浓度的定量反演。误差计算结果说明，该模型的相对误差较低，均方根误差远远低于叶绿素a浓度的平均值。在2009年的验证中，最大相对误差约为58%，远高于2006年的计算结果，这是由于该数据对应的19号点位于影像条带范围内，影像值与实测值对应位置略有偏差。此外，遥感数据与实测数据的获取日期存在13 d的间隔，也使得2009年的验证精度略低于2006年的结果，但总体而言，2个时相的验证精度基本满足目前水质参数的估算要求[26]。综上所述，该模型适用于同一区域不同时相的遥感影像，反演精度较高，普适性较强。

4 巢湖叶绿素a浓度时空变化分析

定量遥感技术在叶绿素a浓度时空分布方面的研究颇具实用性，通过在不同时空尺度下叶绿素a浓度的分布规律，可以更深入地理解水体浮游植物的动态变化过程及其与物理过程的耦合机制[31]。利用本文模型2016年6—11月巢湖叶绿素a浓度反演结果如图6所示。

(a) 2016年6月23日 (b) 2016年7月25日

(c) 2016年8月10日 (d) 2016年11月14日

4.1 空间分布特征

由图6结果可见，巢湖富营养化区域主要分布在巢湖西北部和东部湖域，尤其是夏季时节，巢湖东部双桥河和裕溪河入湖流域会出现较大面积的水华暴发，叶绿素a浓度较低区域主要分布在巢湖中部地区。入湖口的水中悬浮物浓度较高，光合作用受限，于是藻类出现暴发性生长。

4.2 季节分布特征

对比分析可知，6月下旬—7月期间，随着夏季温度升高，藻类进入迅速繁衍期，巢湖整体富营养化程度加剧； 7月底—8月，湖体富营养化达到四季顶峰，中度及以上富营养化程度水体面积高达87%； 8月中旬，藻类开始出现衰退。但约9%轻度富营养化水体区域出现一定程度恶化，这可能是由于不同藻类的潜伏期和繁衍期存在差异； 11月藻类逐渐衰退，大部分水体恢复常态。

由于温度、光照、降水量和营养盐等因素的差异，藻类在水中进行季节性运动。通常藻类于春季休眠，夏季生长暴发，秋季逐渐消亡直至冬季恢复休眠[32]。此外，叶绿素生长规律与模型反演结果大致相似，也在一定程度上证明了本模型在同时反演内陆湖多时相叶绿素a浓度方面的优势及可行性。

5 结论

1)本文基于辐射传输机理，建立叶绿素a浓度反演模型。经巢湖2006年及2009年2组数据验证，得出模型决定系数分别为0.877 8和0.848 8，均方根误差远低于平均值，证明该模型适用于同时反演内陆湖不同时相的叶绿素a浓度，普适性较强。进而应用模型反演并分析了巢湖2016年叶绿素a浓度时空分布特征，与藻类生长规律相一致。

2)以往巢湖叶绿素a浓度反演模型大多依赖实测光谱的统计关系，具有时相局限性。本模型计算过程完全依托影像完成，适用于内陆湖不同时相数据的同步反演。

3)模型普适性的另一表征是不同水域的可移植性。但是，本模型不适用于浅海水域的水质参数反演。后期可以此为研究方向，进一步提高模型的普适性。