基于国产卫星多光谱影像的河流水体浊度遥感联合反演研究

肖 潇，徐 坚,2，3，赵登忠，3，赵保成，3，徐 健，3，程学军，3，李国忠，3

(1.长江科学院 空间信息技术应用研究所，武汉 430010； 2.武汉大学 水利水电学院，武汉 430079； 3.武汉市智慧流域工程技术研究中心，武汉 430010)

1 研究背景

浊度是河流水环境监测重要水质参数之一[1]，且与悬浮物之间存在密切联系[2-6]。对浊度变化分析有助于了解水体总悬浮物或沉积物分布情况，进而为研究水体中污染物沉降、分解和扩散等行为提供有效信息。因此，监测浊度空间分布情况对于水体污染管控具有重要意义。

浊度时空分布特异性很强，具有强季节变化特性[7]，常规采样加实验室分析法不能准确表征浊度时空变化情况。因此，将遥感技术应用于浊度估算，可为研究人员提供其时空动态变化情况。与其他水质参数(如叶绿素a、悬浮物、溶解性有机物)相比，浊度遥感反演研究开展较晚，但也取得了一定成果。基于MODIS数据、TM数据、ERS-2 SAR数据等作为遥感数据源，研究人员在芬兰湖[8]、芬兰湾[9]、法国Adour河口[10]、巴西多斯河流域[11]等区域开展水体浊度遥感反演研究，结果显示反演精度尚可。我国研究人员在浊度遥感反演研究中采用多光谱数据和实测高光谱数据，在我国查干湖[12]、官厅水库[13]、北部湾[14]、珠江口[15]、东昌湖[16]、汉江干流[17]、东中国海[18]、太湖流域[19]、长江口及东海[20]等水域，基于经验法和机器学习法构建了水体浊度反演模型，并根据反演结果分析研究区浊度空间分布情况。

前期研究成果表明，多光谱遥感数据由于具有较高空间分辨率、时间分辨率且更易获取等特性，在实际应用中被广泛采用，成为水体水质遥感业务化监测主要数据来源。在反演方法方面，目前常见的浊度遥感反演研究多通过经验法构建单一模型进行直接反演，但这些单一模型在模拟与预测中表现出的精度不一致，导致各单一模型的反演能力均存在一定局限性、实际应用存在不确定性，反演精度或延展性往往不够理想[21]。因此，集合各模型优势、提高模型精度、降低模型在实际应用中的不确定性具有重要意义。研究表明，起源于气象预报领域的集合建模通过综合各模型信息，可有效提高模型模拟精度和模型稳定性[22-27]，综合利用各种模型所提供的信息而构建的联合反演模型将在模拟精度、模型稳定性及扩展性上较单一模型更具优势。

基于此，本文选择汉江中下游典型河段为研究区，利用2012—2013年原位观测数据，以具有较高时间分辨率和空间分辨率的国产卫星数据作为多光谱遥感数据源，在水体光谱特征和浊度敏感波段分析的基础上，综合利用各种单一模型信息，构建水体浊度多光谱遥感联合反演模型，分析模型反演精度、适用性及稳定性，并基于反演结果分析水质参数浓度时空分布特点。

2 研究区与数据

2.1 研究区

汉江是长江中游最大支流，干流主要位于湖北省境内，全长652 km，流域面积6.4万km2，约占湖北省国土面积的40%[28-30]。其中丹江口水库至钟祥河段为中游，长约270 km，流经低山丘岗，水量和含沙量大增，河道不稳定；下游为钟祥以下河段，长约382 km，流经汉江平原，水流交换，曲流发育，汛期洪水在此聚集，宣泄不畅，易发洪涝灾害。

改革开放以来，汉江流域经济迅速发展，汉江中下游周边如襄阳、钟祥等城市工业发展迅猛，出现了大量污染严重的中小企业，加上固有老污染企业，工业废水排放量不断增加，使中下游河段受到不同程度污染。且伴随南水北调中线工程运行，汉江中下游江段流量减少，流速变缓，水体对沿岸污染物的稀释、自净能力变弱，区域水环境面临着更大挑战。

2.2 数 据

2.2.1 原位观测数据

根据长时间序列汉江中下游水体水质变化情况及野外实地调查结果，预设汉江中下游典型河段采样区域(潜江至仙桃段)、采样点。水质数据采集日期为2012年3月26—27日、5月16—17日、11月19—20日，2013年4月27—28日、6月13—14日、11月18—19日，浊度通过便携式水质分析仪测得。使用美国ASD公司ASD- fieldspec4野外便携式光谱测量仪采集水体高光谱数据；水体光谱测量采用水面以上测量法[31-33]。采样点布设见图1，实测水质数据统计见表1。

表1 实测浊度统计Table 1 Statistics of turbidity in Hanjiang River

图1 汉江研究区采样点分布Fig.1 Locations of sampling sites in Hanjiang River

2.2.2 卫星遥感数据

根据各研究区采样时间及国产卫星数据获取情况，本文以汉江中下游典型河段为研究区，选用HJ-1A/1B CCD数据、GF-1 WFV数据和ZY-3 MUX数据作为遥感数据源(相关传感器基本参数可从中国资源卫星应用中心网站(http:∥www.cresda.com/CN/)获得。本文使用的多光谱影像数据与采样日期对应情况见表2。

表2 多光谱影像数据与采样日期对应情况Table 2 Multi-spectral image data and corresponding sampling time

3 模型与方法

3.1 CM-BP模型

集合建模的关键在于集合方法的选择，即各模型权重确定方法，常用集合方法有平均法[34-38]、熵权法[39-40]、集对分析法[41]、遗传算法[42-43]、神经网络[44-46]和粒子群算法[47-49]等。通过对比常用集合方法可知，神经网络算法特有的非线性映射和自适应学习能力优势，可有效结合各单模型有效信息，既可以弱化单一模型对预测结果的影响，又可以避免预测结果对单一模型的过度依赖，能最大程度发挥各模型优势；同时神经网络算法能够快速完成大量运算，且使用方便、操作简单。

本文从单波段和多波段组合等角度引入多种单一模型，进行单一模型适用性分析，以相关性较大(R>0.6)的单一模型作为联合模型输入。对比常见的集合方法，综合考虑模拟精度高、使用方便、操作简单等因素，选择BP神经网络算法作为集合方法构建内陆水体水质参数多光谱联合模型(CM-BP)。根据联合模型通用式(1)可知，求解最优神经网络结构的过程即为确定最优权重ωi的过程。联合模型通用表达式为

(1)

式中:Y表示联合模型;Xi表示参与联合模型构建的单一模型;ωi表示赋予每个单一水质参数反演模型Xi的权重；N为模型数。

此外，为补充部分同步遥感数据缺失问题，同时也为降低不同传感器设置与大气校正对反演精度的影响，本文根据2013年实测高光谱数据以及被选用的传感器光谱响应函数模拟多光谱数据，再结合野外实测水质参数数据，构建基于多光谱数据的汉江水体浊度联合反演模型。

本文以适当方式筛选、集合多种水质参数反演模型，综合利用各种模型所提供的信息而构建的联合反演模型将在模拟精度、模型稳定性及扩展性上较单一模型更具优势。

3.2 波段匹配分析

本文使用2012—2013年实测水质及光谱数据，采用变量投影重要性指数法提取研究区水体浊度敏感波段，该指数计算公式为

式中:VIP是变量投影重要性指数；p是自变量个数；m是偏最小二乘法从原变量中提取的成分个数；th代表第h个成分；R(Y,th)代表成分th对因变量Y的解释能力，为二者相关系数的平方；whj是轴wh的第j个分量，用于测度第j个自变量对构造th成分的边际贡献。已知第j个自变量对Y的解释通过th传递，如果th对Y的解释能力很强，而第j个自变量在构造th时又起到了重要作用，就认为该自变量对Y的解释能力被视为很大，它在解释因变量Y的时候具有更加重要的作用[50]。

通过筛选，研究区水体浊度敏感波段为549、560～607、615～645、673～677、688～693、772～778、817～819、821～822、824 nm等波段区域。根据上述筛选结果，讨论各多光谱数据波段与浊度敏感波段的匹配情况，分析各卫星数据适宜性。敏感波段与HJ-1A/B CCD数据、ZY-3 MUX数据和GF-1 WFV数据波段匹配情况见图2。

图2 浊度敏感波段与多光谱数据波段匹配情况Fig.2 Matching of sensitive bands and multi-spectral data bands for turbidity

从图2可以看出，HJ-1A/B CCD数据、GF-1 WFV数据和ZY-3 MUX数据波段设置类似，且基本覆盖了浊度的敏感波段。因此，这几种多光谱数据可用于开展浊度反演。根据匹配情况，选用B2、B3、B4波段构建浊度反演模型。

3.3 光谱重采样方法

根据多光谱传感器各波段响应函数，将实测高光谱数据重采样至多光谱数据。多光谱传感器各波段响应函数计算公式为

(3)

式中:Ri为传感器第i波段反射率;φi(λ)为第i波段在波长λ处的响应函数;r(λ)为波长λ处的反射率。

3.4 模型精度验证方法

利用实测水质数据及多光谱数据构建水质参数反演模型，其中3/4样本用于建模，剩余1/4样本用于验证。模拟结果的优劣及模型性能采用决定系数(R2)、均方根误差(RMSE)、相对误差(RE)和平均相对误差(ARE)进行评价。相关计算公式见式(4)—式(7)。

(4)

(5)

(6)

(7)

4 结果与讨论

4.1 CM-BP模型构建与分析

以HJ-1A/B CCD数据、GF-1 WFV数据和ZY-3 MUX数据B2、B3、B4波段反射率模拟值作为反演模型建模参数，从单波段、多波段组合等角度引入多种单一模型。分析这些单一模型与水质参数浓度之间的相关性，挑选相关系数较高(R>0.6)的单一模参与联合反演模型构建，即作为BP神经网络输入层。基于相关性分析，选取相关性>0.6的组合作为BP神经网络输入，对应样本浊度作为输出层神经元。

3种卫星数据的入选波段组合形式分别如下：HJ-1A/B CCD数据入选波段组合形式为B3/B2、B2/B3、(B3+B4)/B2、B2/B4、B3/B4、(B3+B4)·ln(B2)等6种，GF-1 WFV数据入选波段组合形式为B3/B2、B2/B3、(B3+B4)/B2、B2/B4、(B3+B4)·ln(B2)等5种，ZY-3 MUX数据入选波段组合形式为B3/B2、B2/B3、(B3+B4)/B2、B2/B4 这4种。以上述波段组合形式为BP神经网络输入层，经过多次迭代及网络结构调整，求得最优神经网络结构。求解最优神经网络结构，即确定网络层数、神经元数目、每层激励函数和学习速率等BP神经网络基本参数。在进行训练之前，通常会先确定网络层数，再通过调整隐含层中的神经元数目来提高计算精度，这样不仅容易调整，训练效果也更易观察，从BP神经网络构成可知，多层网络适用于解决非线性问题，但过多网络层数又会使网络复杂化，延长网络训练的时间；而在激励函数和输出层函数的选择上，研究人员多采用具有非线性放大系数功能的S型函数为激励函数，线性激活函数为输出层函数。据此，本文使用包含一个隐含层的BP神经网络(即网络层数为3层)，根据网络结构及实际情况，通过调整神经元个数，训练函数和其余参数来求得最优网络结构。为使不同算法拟合结果具有可比性，确定学习速率为0.05，目标误差为0.000 5。

从神经网络训练结果可看出，基于不同卫星传感器的网络的基本参数和训练次数不尽相同，但隐含层神经元个数相近(5个或6个)，训练函数均为traingdx函数，可认为traingdx函数的拟合精度高、适用性强，性能稳定。从拟合结果可看出，基于不同卫星传感器的反演精度均>0.9，平均误差均<10%，拟合精度高；在浊度较低时，反演误差均有增大，但结果仍在可信范围内，表明CM-BP模型对于不同浓度的水质参数适应性较好。BP神经网络具体参数见表3。CM-BP模型拟合情况见图3。

表3 CM-BP模型结构及精度Table 3 Structure and precision of CM-BP model

图3 CM-BP模型浊度拟合情况Fig.3 Fitting of CM-BP model for turbidity

4.2 波段组合模型构建与分析

以相关性最高的波段组合构建浊度反演模型，从拟合结果来看，参与建模的波段组合均为B3/B2，这主要是因为3种卫星传感器的波段设置类似；波段组合模型的拟合精度低于CM-BP模型，且在浊度较低的情况下，拟合误差变大，这是由于波段组合模型是基于统计回归构建的，模型对于参与计算的样本点有较大的依赖，稳定性略差。浊度拟合情况见图4。拟合公式如表4。

图4 基于波段组合模型浊度拟合情况Fig.4 Fitting of combined-band models for turbidity

表4 浊度多波段组合反演模型及精度Table 4 Band-combined models and their precision

4.3 模型验证及评价

利用剩余样本对上述模型进行验证。从建模结果及验证结果可以看出，在反演精度方面，CM-BP模型的精度更高(R2均在0.9以上)，波段组合模型建模精度对应决定系数R2仅在0.6左右，验证精度决定系数在0.8左右；在适应性方面，CM-BP模型优于波段组合模型，主要表现在当浊度浓度较低时，两种模型的反演误差均增大，CM-BP模型的反演结果仍在可信范围内，而波段组合模型的反演结果不可接受；在模型稳定性方面，CM-BP模型的建模和反演精度较为一致，波段组合模型的预测精度高于建模精度，说明模型对于参与计算的样本点有较大的依赖，稳定性较差。综合比较两种模型，认为CM-BP模型性能更优。验证结果见表5。

表5 模型验证结果Table 5 Verification of retrieval models

4.4 汉江中下游浊度时空分布

根据前述构建的浊度联合反演模型，使用同步多光谱影像开展研究区水体浊度模拟及时空变化分析。根据所构建的浊度联合反演模型，使用多源同步多光谱数据开展研究区水体浊度模拟，并制作浊度时空分布专题图(见图5)。反演结果表明，春季采样期，研究区水体浊度主要处于20.00～50.00 NTU范围，潜江河段浊度略低于仙桃河段，主要处于20.00～40.00 NTU范围，仙桃河段主要处于30.00～50.00 NTU范围；夏季采样期，研究区浊度略高于春季，水体浊度主要处于40.00～60.00 NTU范围，其中潜江河段略高于仙桃河段，主要在45.00～60.00 NTU之间，仙桃河段在45.00 NTU左右；秋季采样期，水体浊度主要处于15.00～25.00 NTU范围，显著低于春、夏两季，潜江河段浊度略高于仙桃河段。对比反演结果与实测数据，两者较为符合，说明CM-BP可以较为准确地反映浊度的空间分布状况。

图5 浊度时空分布情况Fig.5 Temporal-spatial distribution of turbidity

此外，从图5 (b)—图5(f)可以看出，基于不同时空分辨率多光谱影像数据得到的浊度的空间分布在整体上趋势基本一致，说明CM-BP联合反演模型可同时适用于3种不同的卫星数据，对不同时空分辨率数据的适应性较强。基于3种卫星遥感数据的反演结果虽在整体趋势上基本一致，但从反演结果可以看出(以仙桃段为例，“W”型河段)，基于HJ-1A/B CCD数据在河流水体边界处的反演结果偏低，而基于ZY-3 MUX和GF-1 WFV数据未出现这种情况。这主要是由于所选取的研究区域河段宽度最窄处不足150 m(平水期)，水体宽度较窄，HJ-1A/B CCD数据较低的空间分辨率将影响河流水体边界的准确提取，进而影响水体浊度反演精度，导致最终边界处反演结果偏低。空间分辨率的提高在一定程度上可以提高反演精度，故空间分辨率亦是一个重要的数据源选择因素。

5 结论与展望

本文引入集合建模思想，通过相关分析和神经网络算法对单一模型进行分析、联合，构建基于卫星多光谱影像数据的内陆水体水质参数联合反演模型(CM-BP)，以汉江中下游典型河段为研究区，利用实测水质、高光谱数据，以及国产卫星多光谱数据开展浊度遥感反演研究，得到如下主要结论：

(1)CM-BP模型反演精度显著高于波段组合模型反演精度；CM-BP模型适应性更强，当浊度浓度较低时，浊度反演误差增大，但CM-BP模型反演结果仍在可信范围内，而波段组合模型反演结果不可接受；模型稳定性方面，CM-BP模型建模和反演精度较为一致，波段组合模型预测精度和建模精度存在一定差别，说明波段组合模型对于参与计算样本点有较大依赖，稳定性较差。CM-BP模型基于浊度的敏感波段和非线性算法构建，参与建模波段不受人为主观影响，比常规经验算法选取的波段更具可信度，提高了建模效率。综合考虑模型反演精度、适应性、稳定性和可信度，认为CM-BP模型性能更优。

(2)国产卫星如环境与灾害监测预报小卫星星座系统、高分一号卫星、资源三号卫星的多光谱数据，可以满足高精度、实时及大尺度的内陆水体水质遥感监测工作需求，均可作为河流水体水质遥感反演研究及水质遥感监测业务化工作的优选数据源。而从光谱分辨率、时空分辨率以及影像幅宽等因素考虑，高分一号卫星WFV数据在河流水体水质遥感反演研究上更具优势。

(3)遥感数据空间分辨率对反演精度与误差的影响，尤其是在水陆交界区域，主要与河流水体范围是否被精确识别有关。因此，在遥感数据空间分辨率满足水体精确识别的条件下，只要所选择的遥感影像数据分辨率和波段设置适用于水质反演，本文构建模型均可适用。

(4)影响河流浊度的因素诸多，如含沙量、突发性水污染、流速等，厘清这些因素对反演精度的影响是未来进一步研究的关键内容。受采样次数的限制，本文针对上述因素对反演结果影响的讨论尚不够深入，后续将继续积累数据进行进一步研究。