采用芬兰的机载高光谱成像系统(airborne hyperspectral imaging systems, AISA)获取了研究区域的航空高光谱数据, AISA获取数据的光谱范围为400~970 nm, 波段数126, 航高1 000 m, 地面空间分辨率约为1.2 m, 航测区域范围约80 km², 共获取航线14条。 水面和陆面测量实验是地面同步实验的主体, 采样点分布根据水体分布情况确定, 共采集75个样点, 如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区域采样点分布图Fig.1 Distribution of sample points in study region

2012年11月5日— 7日和11月19日分4架次获取航空高光谱数据。 2012年10月26日, 11月1日, 11月5日— 6日期间, 在晴空无云、 水体较稳定的情况下选取10至15时时间区段, 进行地面数据的采集。 采集的数据包括: 水体样本、 水面以上光谱数据及采样空间位置等辅助数据。 水面光谱测量使用ASD 便携式野外光谱仪FieldSpec® Pro FR(波长范围350~2 500 nm), 采用“ 水面以上法” 测量水面光谱, 测量时仪器观测平面与太阳入射平面夹角135° , 探头与水面法线夹角40° 的观测角度。 为保障地面数据的代表性, 选取了不同类别区域进行水质采样, 包括湖泊、 湿地、 农业用地河网、 航运河网和城市内部河网的水体进行采样; 同时选择开阔的陆面区域, 测量典型地物: 植被、 裸土、 道路等光谱数据用于航空高光谱数据的大气校正, 水面采集的水样送到水质检测站进行分析处理, 叶绿素a采用分光光度法测定, 悬浮物采用重量法测定, 获得2项水质参数浓度实测数据。

航空高光谱数据预处理包括辐射校正、 几何校正、 大气校正、 航带拼接以及水体提取, 其中对水质参数反演最重要的是大气校正。 大气校正实现图像像元从表观辐亮度到地表反射率的转换, 大气条件输入的准确性直接影响地表反射率反演精度, 采取辐射传输模型法进行大气校正, 根据具体的成像时间和地理位置, 大气模式选择中纬度冬季模式, 气溶胶模型选择城市型, 能见度和水汽含量从附近气象站数据获取, 将这些大气条件参数等输入大气辐射传输模拟软件MODTRAN中计算, 获得地表反射率数据。 进一步, 对研究区域水体提取实现水陆分离, 获取嘉兴的河网水体范围。

通过对水面同步实测的高光谱数据分析, 利用短波红外处水表反射率数据, 考虑水体在短波红外反射率应该为0, 将不符合的异常数据剔除, 最终确定有效数据为50个, 利用式(1), 获取离水辐亮度, 最终获取水表反射率光谱R_rs曲线, 如图2所示。

Lw=Lsw-rLsky(1)

式中, L_w为离水辐亮度, L_sw为测量获得的水面以上辐亮度, r为气-水界面反射率, L_sky是不包含水体信息的天空光。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 实测水体反射率光谱Fig.2 The measured water reflectance spectra

由水表反射率光谱曲线可见, 嘉兴市城市河网水体的水表反射率光谱曲线整体有4个比较明显的反射峰, 分别为580 nm附近、 650 nm附近、 700 nm附近和820 nm附近, 其中前三个反射峰组成了一个大的反射峰, 水表反射率从580~700 nm呈下降趋势, 600~650 nm水表反射率变化较小, 700 nm附近光谱反射率开始急剧下降, 850 nm之后由于纯水的吸收作用, 水表反射率呈现趋于0的变化。

1.3.1 理论基础

水体辐射传输模型是水体光谱特性分析的理论基础, 水体辐射传输模型是进行水体光学遥感建模的重要工具, 它模拟了水体辐射传输过程。 水色遥感中一般不直接解算水体辐射传输方程, 而是经常使用由辐射传输方程经过简化近似后得到的生物光学模型。

实现水体固有光学量的反演首先需要构建表观光学量与固有光学量之间的联系, 模拟水体中的辐射传输过程, 相关学者做了大量的研究, 给出了吸收系数、 后向散射系数与水表反射率的关系模型。

1.3.2 QAA算法和IIMIW算法实验研究

目前常用于内陆Ⅱ 类水体固有光学量反演算法有QAA算法、 IIMIW算法等。

QAA算法是由Lee等^{[1, 2]}于2002年提出并不断发展的一种反演水体固有光学量的算法, QAA方法是基于水表反射率和IOPs半经验关系的有效方法, 并且已经被近海岸和海洋区域获取的数据验证。 QAA算法的处理过程主要分为两个部分, 第一部分完成吸收系数和后向散射系数的反演, 第二部分利用反演获得的总吸收系数反演水体主要组成成分的吸收系数。 利用QAA算法对实测水面以上光谱进行固有光学量反演, 获得剔除纯水总吸收系数, 如图3(a)所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 剔除纯水的总吸收系数 (a): QAA算法; (b): IIMIW算法Fig.3 Total absorption coefficient eliminated pure water (a): QAA algorithm; (b): IIMIW algorithm

IIMIW算法是由Li等^[6]提出的一种应用于内陆水体固有光学量的反演算法, 该算法考虑到目前算法在构建生物光学模型中的常量应该为比较复杂的变量, 然后构建了波段比值的方式进行了消除, 同时利用经验模型获取参考波段的后向散射系数, 最后计算全部波段的散射系数和吸收系数。 利用IIMIW算法对实测水面以上光谱进行固有光学量反演, 如图3(b)所示。

1.3.3 改进QAA算法

分析应用QAA和IIMIW算法计算内陆水体存在一定的问题, 首先两种算法使用的Y系数估计模型都是经验模型, Y系数是利用NOMAD数据集(NASA bio-optical marine algorithm data set)计算得出, 这些数据都是开阔海洋或近海水体, 与实际内陆水体有较大的差别, 因此这些经验模型并不一定适用。 第二, QAA算法是基于海洋遥感发展的, 含有大量悬浮物和有色溶解有机物内陆浑浊水体光谱更为复杂, QAA算法不一定适用; 再则, QAA算法在估计总吸收系数的参考波段时是通过模拟数据得出, 所选取的参考波段λ ₀为555 nm不适用于内陆水体, IIMIW估计后向散射系数的参考波段为560 nm, 与QAA算法的参考波段十分接近。 最后, 对于IIMIW算法利用波段比值最终计算其他波段的吸收系数时选取了709 nm, 假设709 nm处的吸收系数只有纯水作用, 最终反演得到的参考波段后面的数据除了纯水吸收系数以外远大于0, 出现了数据异常, 而且河网水体可能存在泥沙等影响, 在红光近红外水表反射率增加, 这个假设在河网水体基本不成立, 所以直接应用后向散射估计时, 在参考波段后波段吸收系数会出现较大的异常。

通过以上分析并结合河网城市水体特点, 对QAA算法做了适应性改进。 通过表1和表2对比, 可以看出步骤5的后向散射估计模型进行了调整, 步骤2— 步骤5的参考波段转向近红外波段。

表1

Table 1

表1(Table 1)

表1 Lee QAA算法具体步骤 Table 1 Steps of Lee quasi-analytical algorithm 步骤 变量 数学公式 方法 0 rrs =Rrs/(0.52+1.7Rrs) 半分析法 1 μ (λ ) = -g0+[(g0)2+4g1rrs(λ)]1/22g1; g0=0.089 5, g1=0.125 半分析法 2 a(λ 0) =aw(λ 0)+1 0-1.146-1.366χ-0.469χ2, χ =log rrs(443)+rrs(490)rrs(λ0)+5rrs(667)rrs(490)rrs(667), λ 0=555 经验法 3 bbp(λ 0) = μ(λ0)a(λ0)1-μ(λ0)-bbw(λ 0), λ 0=555 分析法 4 Y =2.0 1-1.2exp-0.9rrs(443)rrs(λ0), λ 0=555 经验法 5 bbp(λ ) =bbp(λ 0) λ0λY, λ 0=555 半分析法 6 a(λ ) = 1-μ(λ)bbw(λ)+bbp(λ)μ(λ) 分析法

表1 Lee QAA算法具体步骤 Table 1 Steps of Lee quasi-analytical algorithm

表2

Table 2

表2(Table 2)

表2 改进QAA算法具体步骤 Table 2 Steps of advanced quasi-analytical algorithm 步骤 变量 数学公式 方法 0 rrs =Rrs/(0.52+1.7Rrs) 半分析法 1 μ (λ ) = -g0+(g0)2+4g1rrs(λ)1/22g1; g0=0.089 5, g1=0.125 半分析法 2 a(λ 0) =aw(λ 0), λ 0=812 经验法 3 bbp(λ 0) = μ(λ0)a(λ0)1-μ(λ0)-bbw(λ 0), λ 0=812 分析法 4 Y =2.0 1-1.2exp-0.9rrs(443)rrs(λ0), λ 0=812 经验法 5 bbp(λ ) =bbp(λ 0) λ0λY+ρ Rrs(λ 0), ρ =10, λ 0=812 半分析法 6 a(λ ) = [1-μ(λ)][bbw(λ)+bbp(λ)]μ(λ) 分析法 注: Rrs为水面以上反射率; rrs为水面以下反射率; μ 为水面以上反射率与水体总吸收系数和后向散射系数的关系变量; g0, g1和ρ 为经验系数; λ 为波长; λ 0为参考波段; a为水体总吸收系数; bbp为颗粒物后向散射系数; bbw为纯水后向散射系数; Y为颗粒物后向散射系数随波长的衰减指数 Note: Rrs is the remote sensing reflectnace above water surface; rrs is the remote sensing reflectance below water surface; μ is a variable to describe the relationship between the remote sensing reflectance below water surface with total absorption and backscattering coefficients of water column; g0 and g1 and ρ are empirical coefficients; λ is wavelength; λ 0 is the reference wavelength band; a is the total absorption coefficient of water column; bbp is the backscatter coefficient of the suspended particles; bbw is the backscatter coefficient of the pure water; Y is a attenuation index to characterize the attenuation of backscatter coefficient with the wavelength for the particles

表2 改进QAA算法具体步骤 Table 2 Steps of advanced quasi-analytical algorithm

步骤5颗粒物后向散射估计模型调整。 由于河网水体相比其他常规水体(海洋、 湖泊等), 水体流动性导致水体中的泥沙等悬浮物浓度较高, 尤其在高悬浮物含量的水体中后向散射系数并不满足幂指函数的关系, 河网水体相比其他水体的后向散射多出的散射由无选择性散射主导, 那么水体的后向散射受波长影响较小, 在假设泥沙等悬浮物浓度对光谱反射率随波段变化较小情况下, 可以得出式(2)经验估计模型。

bbp(λ)=bbp(λ0)λ0λY+ρRrs(λ0)(2)

式中, b_bp(λ )和b_bp(λ ₀)分别为所有波段和参考波段颗粒物后向散射系数, ρ =10为经验系数, Y为颗粒物后向散射系数随波长的衰减指数, 为参考波段的水表反射率。

步骤2— 步骤5参考波段调整。 已有的研究证实应用于浑浊水体时参考波段的选择应转向更长的谱段, 如果使用NIR(near infrared)波段, 总吸收系数可以基本被纯水的吸收系数替代, 由于水的吸收系数远远大于其他物质, a(λ ₀)≈ a_w(λ ₀), 其中a为总吸收系数, a_w为纯水吸收系数。

使用改进的QAA算法对实测水面以上光谱进行固有光学量反演, 得到浮游植物吸收系数[图4(a)]和剔除纯水总吸收系数[图4(b)]。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 改进QAA算法反演水体固有光学量 (a): 浮游植物吸收系数; (b): 剔除纯水总吸收系数Fig.4 Inherent optical quantity retrieved by improved QAA algorithm (a): Absorption coefficient of phytoplankton; (b): Absorption coefficient eliminated pure water

对于QAA算法获取剔除纯水吸收系数从图3(a)可以看出, 吸收系数整体偏小, 吸收系数在700 nm之前呈现下降趋势, 在700 nm后的数据基本是异常值, 而且在650~700 nm之间的吸收系数出现了负值, 表明实际反演的总吸收系数也较低。 而对于图4(b)中改进QAA算法反演的剔除纯水吸收系数结果, 考虑无选择性散射和参考波段为812 nm时, 在参考波段前面的负值异常值以及近红外吸收系数较大的异常值基本消除, 具有较好的效果, 符合水体的吸收系数实际情况。

浮游植物通常指水中以浮游状态生活的小植物, 叶绿素a含量与浮游植物的含量等信息直接相关。 悬浮物通常指不溶于水的泥沙、 粘土、 藻类等, 悬浮物含量与剔除纯水后的含量直接相关。 计算叶绿素a浓度与浮游植物吸收系数波段组合, 以及悬浮物浓度与剔除纯水吸收系数波段组合的相关系数, 获得最佳波段组合, 如图6所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 水质浓度与吸收系数波段组合的相关系数 (a): 叶绿素a浓度和浮游植物吸收函数; (b): 悬浮物浓度和剔除纯水吸收系数Fig.5 Correlation coefficient of concentration of water quality & combination of absorption coefficient band (a): Concentration of Chlorophyll⁃a & absorption coefficient of phytoplankton; (b): Concentration of suspended solids & absorption coefficient eliminated pure water

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 吸收系数最优波段比值与水质浓度线性回归分析 (a): 浮游植物吸收系数和叶绿素a浓度; (b): 剔除纯水吸收系数和悬浮物浓度Fig.6 Linear regression analysis of optimal band ratio between absorption coefficient & concentration of quality parameters (a): Absorption coefficient of phytoplankton & concentration of Chlorophyll-a; (b): Absorption coefficient of eliminated pure water & concentration of Suspended Solids

图5(a)中叶绿素a浓度与浮游植物吸收系数的双波段比值相关系数位于0.63以上, 呈现显著相关, 而且叶绿素a浓度与浮游植物吸收系数波段组合的相关系数既有正值也有负值。 叶绿素a浓度与中心波长为677.2和725.1 nm的浮游植物吸收系数的相关性最好, 因此选择最优波段为678和726 nm, 这两个波段比值与叶绿素a浓度最大相关系数为0.75, 具有很强的线性相关性。

图5(b)中悬浮物浓度与剔除纯水吸收系数的双波段比值相关系数位于0.65以上, 呈现显著相关, 而且悬浮物与剔除纯水吸收系数波段组合的相关系数既有正值也有负值。 悬浮物浓度与中心波长为734.6和793.0 nm的剔除纯水吸收系数的相关性最好, 因此选择最优波段为734和790 nm, 这两个波段比值与悬浮物浓度最大相关系数为0.82, 具有很强的线性相关性。

浮游植物吸收系数最优波段比值B726/B678与叶绿素a浓度线性回归分析, 剔除纯水吸收系数最优波段比值B734/B790与悬浮物浓度线性回归分析, 如图6所示。

图6(a)得到叶绿素a浓度反演模型见式(3), 拟合结果的决定系数R²为0.64, 具有较强的线性相关性; 均方根误差RMSE为1.76。

c(Chl⁃a)=8.717×aph(726)aph(678)-3.068(3)

式中, c(Chl-a)为叶绿素a浓度, a_ph为浮游植物吸收系数。

图7(b)得到悬浮物浓度反演模型见式(4), 拟合结果的决定系数R²为0.71, 具有较强的线性相关性; 均方根误差RMSE为15.63。

c(SS)=66.05×at⁃w(734)at⁃w(790)-62.8(4)

式中, c(SS)为悬浮物浓度, a_t-w为剔除纯水吸收系数。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 高光谱数据反演水质参数浓度分布图 (a): 叶绿素a浓度分布; (b): 悬浮物浓度分布Fig.7 Concentration distribution of water quality parameters retrieved using hyperspectral date (a): Distribution of Chlorophyll-aconcentration; (b): Distribution of suspended solid sconcentration

利用研究区航空高光谱数据和上述水质参数反演模型, 反演得到研究区叶绿素a和悬浮物的浓度分布, 如图7所示。

针对叶绿素a水质参数高光谱定量反演结果, 利用与航空高光谱数据同步采集且未用于建模的4个水质采样化验数据进行对比分析, 结果如表3所示, 4个样点的平均相对误差为9.2%。

表3

Table 3

表3(Table 3)

表3 叶绿素a水质参数反演和实测对比 Table 3 Comparison of retrieved and measured water quality parameters of Chlorophyll-a 样点 实测值/ (mg· m-3) 反演值/ (mg· m-3) 72 0.480 0.430 73 4.660 5.400 74 35.900 39.100 75 21.600 23.500 平均相对误差/% 9.2

表3 叶绿素a水质参数反演和实测对比 Table 3 Comparison of retrieved and measured water quality parameters of Chlorophyll-a

图7(a)中叶绿素a浓度较高的地方集中在北郊河北部的鱼塘以及西半部分农田中的池塘, 同时在东外环河部分河的段也出现较高的值, 叶绿素a浓度较低的区域集中在城市中央的南湖、 西南湖、 环城河区域, 而且在南湖等大面积水域也表现了叶绿素a浓度由湖边向湖心逐渐降低的趋势, 外围河道区域的叶绿素a浓度略高于城市内部区域, 这与在野外数据采集看到的航道有大量水葫芦等水生植物也比较相符, 在石臼漾湿地区域净水沉淀池叶绿素a浓度也表现了由边缘向中心变小的趋势, 而且湿地出水口的叶绿素a浓度也略低于入水口。

针对悬浮物水质参数高光谱定量反演实验结果, 利用与航空高光谱数据同步采集且未用于建模的4个水质采样化验数据进行对比分析, 结果如表4所示, 4个样点的平均相对误差为9.4%。

表4

Table 4

表4(Table 4)

表4 悬浮物水质参数反演和实测对比 Table 4 Comparison of retrieved and measured water quality parameter of suspended solids 样点 实测值/ (mg· L-1) 反演值/ (mg· L-1) 72 6.000 4.670 73 21.000 25.220 74 42.000 51.740 75 27.000 23.440 平均相对误差/% 9.4

表4 悬浮物水质参数反演和实测对比 Table 4 Comparison of retrieved and measured water quality parameter of suspended solids

图7(b)中悬浮物线性反演结果分布基本呈现以南湖为中心, 向四周扩散浓度逐渐偏高的趋势, 由于外围主航道受到水流速以及船只扰动的结果, 北郊河、 东外环河泥沙含量较高, 悬浮物浓度最高, 城市中部各水系石臼漾湿地、 新塍塘水系西段、 京杭大运河(中部段)、 穆河溪、 长纤塘、 六里长泾水系、 嘉善塘东段悬浮物浓度次高; 南湖、 西南湖、 湘家荡、 环城河、 京杭大运河(西段)水系浓度最低。

根据航空高光谱数据反演的水质参数与地面水样水质浓度对比, 叶绿素a水质浓度反演结果相比悬浮物较好, 与水质反演模型构建时的均方根误差对比一致。 用于分析反演结果的4个地面样点数据, 处于航空高光谱数据同步获取区域, 实验结果分析的条件较好。 4个样点实测值对应的反演值, 是对验证样点位置半径5米反演结果平均值作为对比的反演值, 提高了结果分析准确性。 从叶绿素a和悬浮物浓度分布特征分析上, 也基本与嘉兴市城市河网实际运行情况相符。