基于多时相主被动遥感的漓江水面监测与水质参数反演(2016-2020年)*

刘 曼，付波霖，何宏昌，解淑毓，仇继圣，孙习东，龚烨云，劳植楠，左萍萍

(桂林理工大学测绘地理信息学院，桂林 541006)

河流水文信息动态监测是分析流域水体演变规律的基础，而水质信息的采集则是评价流域水体质量优劣程度、衡量水体初级生产力的前提，及时准确地监测、获取河流水文和水质信息对河流水资源开发管理、水环境保护、水生态修复以及防洪减灾等方面具有重要意义[1-2]. 漓江是桂林风景的灵魂，是中国唯一一条入选的全球最美河流[3]. 然而，近年来漓江水域水量暴涨、干涸等现象突出，影响了桂林人民的正常生活. 快速、准确地掌握漓江水体信息，并依据流域时空演变特点研究相应对策，对漓江流域进行合理开发、规划具有重大意义.

目前，遥感技术已被用于河流水文、水质信息的监测，相较于传统的地面监测方法，其在宏观性、时效性、便利性上具有较大优势. Yang等[4]和王大钊等[5]使用归一化差异水体指数(NDWI)、改进归一化水体指数(MNDWI)与自动提取水体指数(AWEI)等方法，基于Sentinel-2和Landsat 8影像提取了城市地表水分布信息，结果表明Sentinel-2 影像提取的水体信息更为准确，Kappa系数达0.88以上；段秋亚等[6]基于高分一号(GF-1)数据，采用NDWI阈值法，提取了 2 块不同尺度和不同复杂度的代表性水域，与人工解译的水体信息相比较，阈值法提取精度分别达到 97.59% 与99.15%；Yao等[7]利用城市水自动提取方法和资源三号(ZY-3)影像提取了青岛多个地区的地表水信息，与NDWI指数提取结果相比，该方法Kappa系数提高了7%～32%. 因光学影像易受到天气的影响，部分学者则利用微波遥感全天时、全天候监测以及后向散射系数对水体较敏感等特点来进行水体提取研究. Hsu等[8]基于台湾地区四景Sentinel-1雷达影像，利用后向散射系数提取了河流水面宽度，与无人机影像相比，利用雷达影像测量的结果易受阴影的干扰；Liang等[9]基于Sentinel-1合成孔径雷达图像，采用局部阈值化方法提取了Louisiana部分洪水淹没地区的水体轮廓，与传统阈值法相比，局部阈值法在区别水体与非水体方面精度更高，用户精度与生产者精度提高了4%～13%；谷鑫志等[10]提出了水体自动化提取算法，利用两景不同成像模式的GF-3 SAR影像对湖南省东北部水体进行信息提取，结果表明该算法提取精度达到了85%以上，精度优于基于Landsat 8与GF-1影像的NDWI、MNDWI、AWEI指数. 叶绿素浓度遥感反演研究对水体的管理起着积极作用，为验证遥感影像对流域水质信息的监测效果，Ansper等[11]和Toming等[12]利用Sentinel-2数据反演了Estonian湖泊叶绿素a(Chl.a)浓度，证明了Sentinel数据对水体中叶绿素的估算和湖泊时空动态的跟踪的适用性；Nazeer等[13]为监测香港沿海复杂水域的Chl.a浓度，使用Landsat TM 、ETM+影像中蓝、绿、红和近红外4个波段数据及实测叶绿素数据建立了回归模型，相关系数达0.89；Filipponi[14]利用C2RCC算法从Sentinel-2数据中提取意大利波河水体信息，并分析了水体中Chl.a和总悬浮物质(TSM)浓度变化，证明了利用卫星多光谱数据监测河流水质变化的可行性. 综上所述，以上研究主要集中在比较不同遥感水体指数提取水体的精度，少有学者利用长时间序列的主被动遥感影像提取水体信息和反演水体Chl.a、TSM浓度等水质参数，定量分析河流水文和水质信息的年内动态变化规律.

因此，本文以漓江流域为研究区，以2016-2020年Landsat 8 OLI、GF-1和Sentinel-2A多光谱数据和2018年逐月Sentinel-1A IW模式的SAR影像为主要数据源，利用构建的多种主被动遥感水体指数逐月提取漓江水面信息，反演漓江逐月的Chl.a和TSM浓度，构建了漓江流域的278个基本评价单元，并在此单元内基于河岸线发育率、水面变化区域差异值和水体信息变化动态度定量分析漓江水面、Chl.a和TSM浓度的年内、枯水期和汛期、上、中和下游的时空演化特征，为漓江的水环境保护、流域合理开发与保护提供技术支撑和理论依据.

1 研究区概况及数据源

1.1 研究区概况

漓江流域(24°18′～25°41′N，109°45′～110°40′E)位于广西壮族自治区桂林市境内，漓江发源于兴安县华江乡猫儿山老山界东南侧，自北向南流经兴安县、灵川县、桂林市区、阳朔县、荔浦县、恭城县、平乐县等地，全长164 km，属珠江水系西江的支流[15]. 桂林位于低纬度区域，属亚热带湿润季风气候，雨季集中在夏季期间. 漓江年径流量大，但分布不均匀，水位全年变化显著，3-8月为汛期，9月至次年2月为枯水期，且两期差异显著[16]. 枯水期河床外露，严重时船只无法通行；汛期水量暴涨，可淹没两岸河床和滩涂. 桂林山水甲天下，而漓江则是桂林山水的灵魂所在，漓江部分河段有着喀斯特地貌，是广西东北部喀斯特地貌发育的最典型地区，象鼻山、冠岩、黄布倒影等地较为出名，1982年国务院将漓江流域列入“首批国家重点风景名胜区”之一[17]，桂林旅游业得以快速发展. 近年来漓江水量暴涨、干涸现象突出，因受水位涨幅的剧烈影响，桂林旅游业呈现明显的淡旺季，桂林人民的正常生活也受到严重影响，汛期时因洪水造成的经济损失尤为严重. 漓江流域位置如图1所示.

图1 漓江流域位置Fig.1 Location of Lijiang River Basin

1.2 数据源及数据预处理

为探究漓江水体年内动态变化特征，本文获取了漓江流域全年逐月的主被动遥感影像. 光学遥感影像主要包括：①2016-2020年12个月的Sentinel-2A多光谱数据，下载网址为https://scihub.copernicus.eu/；②2016-2018年11个月Landsat 8 OLI和GF-1 WFV多光谱数据，其中Landsat 8 OLI数据是从美国地质调查局USGS网站(https://earthexplorer.usgs.gov/)下载获取，因部分月份出现Landsat 8 OLI影像研究区覆盖不全、云覆盖率较高等问题，本文将从中国资源卫星应用中心网站(http://www.cresda.com)下载的对应月份的GF-1 WFV多光谱数据作为补充数据，而在2016-2018年研究区3月的Landsat 8影像和GF-1影像均被较多云遮挡，无法进行漓江水体的提取. 本文利用Sen2Cor v2.8插件对Sentinel-2A影像进行辐射校正、大气校正处理，再用SNAP软件将多光谱影像重采样为15 m的空间分辨率，并导入到ENVI 5.4软件中对数据进行波段合成、影像拼接和裁剪，生成研究区影像；而基于Sentinel-2A影像对Chl.a及TSM浓度的反演，则直接使用SNAP软件对原数据重采样，并进行C2RCC大气校正，将Salinity参数设置为0.0001；运用ENVI 5.4软件分别对Landsat 8 OLI多光谱影像进行辐射定标和FLAASH大气校正，对全色波段进行辐射校正，并利用Gram-Schmidt Spectral Sharpening算法将校正后的全色与多光谱进行图像融合，生成空间分辨率为15 m的多光谱影像，最后裁剪得到漓江流域影像数据；同时，利用ENVI 5.4软件对GF-1 WFV影像进行辐射定标、大气校正、几何校正、影像拼接与裁剪.

本文在欧空局遥感数据中心网站(http://www.cnblogs.com/)下载了漓江流域2018年12个月IW模式SLC格式的Sentinel-1A影像. 运用SNAP 软件处理Sentinel-1A 影像，先后经过热噪声去除、辐射定标、Burst去除、Refined Lee滤波处理、多视处理、斜距影像转换为地距影像、地形纠正、地理编码和影像裁剪，导出空间分辨率为15 m为ENVI格式的后向散射系数数据. 为了减少相干噪声对准确提取和解译研究区域水体的影响，本文分别将Refined Lee滤波窗口设置为7×7，多视处理视数为4. 基于 30 m SRTM DEM作为参考DEM 对Sentinel-1A影像进行地理编码，生成投影坐标系为 WGS_1984_UTM_Zone_49N. 研究区的主被动遥感数据集以Sentinel-2A多光谱影像为基准影像，其他数据在ENVI 5.4 软件中进行相对配准. 数据源具体信息如附表Ⅰ所示.

水体Chl.a和TSM浓度实测数据来自于广西壮族自治区桂林市水文水资源局提供的2018年逐月数据，用于验证基于C2RCC算法和Sentinel-2A影像反演的逐月Chl.a和TSM浓度. 漓江水面宽度数据验证数据基于Google Earth中高空间分率遥感影像进行解译和矢量化提取. 其中，在解译提取漓江水面宽度数据时，将Google Earth中影像调整为与本文所用遥感影像对应的月份，并从漓江上、中和下游各选取5个漓江水面横截面计算出平均值作为最终漓江水面宽度的验证数据.

2 研究方法

为系统研究漓江水文和水质信息年内动态变化特征，本文基于水利部《河流健康评估指标、标准与方法 V1.0》，以500 m 为步长等分主航道中心线且生成覆盖整个研究区的 278个矩形缓冲区，格网化漓江水体建立278个基本评价单元. 基于基本评价单元分别计算河岸线发育系数(SDI)、水面变化区域差异值(WDr)及水体信息变化动态度(K)3种指标参数，来定量分析漓江水面、Chl.a和TSM浓度的年内、枯水期和汛期上、中和下游的变化规律. 本文采用阈值分割算法对漓江遥感水体指数进行分割运算，提取逐月的水体信息，分别计算了4类光学遥感水体指数和雷达后向散射系数：Landsat 8 OLI和GF-1多光谱影像分别计算了归一化水体指数(NDWI)、改进型归一化水体指数(MNDWI)、增强型水体指数(EWI)，Sentinel-2A多光谱影像计算了归一化差值池指数(NDPI)，Sentinel-1A SAR影像计算了后向散射系数(S). 综合利用多种光学遥感水体指数和雷达后向散射系数进一步构建了主被动遥感加权指数(JQ)，增强水体信息的同时抑制其他非水体信息. 本文采用二类水体区域性近岸海域水色算法(C2RCC)、最大叶绿素指数(MCI)、双波段比值法(DoubleR)及Chl.a反射峰强度(ρchl)4种方法，反演计算漓江水体Chl.a和TSM浓度，进而探究其年内动态变化规律，具体技术流程如图2所示.

图2 研究方法和技术路线Fig.2 Research method and technical route

2.1 基于主被动遥感的漓江水体提取

由于NDWI水体指数对植被与水体有较好的区分能力，MNDWI水体指数适用于增强和提取背景为建成区的水域信息，EWI水体指数则结合了以上两种指数的优点，较好地对植被、土壤、建筑区进行反射抑制，增强了水体亮度，NDPI水体指数不仅可以提取小型池塘和低于0.01 hm2的水体，还可以区分水体与周围环境中的植被信息[18-21]. 因此，本文分别基于NDWI、MNDWI、EWI、NDPI指数和后向散射系数(S)数据进行灰度分割，确定佳阈值分割提取漓江水体信息. 具体计算公式为：

归一化水体指数(NDWI)：

(1)

改进型归一化水体指数(MNDWI)：

(2)

增强型水体指数(EWI)：

(3)

归一化差值池指数(NDPI)：

(4)

以上各式计算结果均浮点型数据，式中，BGreen为多光谱中的绿波波段，BNIR为近红外波段，BSWIR表示短波红外波段.

为了提高水体提取精度，融合了3种指数和后向散射系数S的水体提取优势，本文在对比NDWI、MNDWI、EWI指数与后向散射系数S提取漓江水体效果的基础上，根据水体提取精度确定加权系数，并利用ENVI 5.4软件中波段运算工具进行加权求和运算，得出主被动遥感加权水体指数(JQ)影像，再对其灰度进行分割、确定最佳阈值提取漓江流域逐月水体信息，具体计算公式为：

JQ=a·BNDWI+b·BMNDWI+c·BEWI+d·BS

(5)

表1 主被动遥感加权水体指数中各组成 部分逐月权重系数

式中，a、b、c、d为权重系数，依据试错法进行反复试验所设定，首先通过阈值分割确定NDWI、MNDWI、EWI与后向散射系数4种指数逐月水面提取结果的最适阈值范围，再依据不同月份各类指数分别对漓江水体的提取效果进行权重设定，逐月权重系数值如表1所示，BNDWI为研究区NDWI指数数据，BMNDWI为研究区MNDWI指数数据，BEWI为研究区的EWI指数数据，BS为研究区后向散射系数数据.

Chl.a是水生植物的重要组成成分，悬浮颗粒物是水体中营养物质和污染物的重要载体，两者皆是反映水体质量优劣程度的重要参数. 基于神经网络技术的C2RCC算法可对多种数据进行大气校正，该算法可输出校正后的遥感反射率、固有光学量、Chl.a浓度、TSM浓度和黄色物质等多种产品数据，本文则基于此算法反演漓江水体Chl.a及TSM浓度[22-23]；Sentinel-2A中MSI传感器有3个波段可获取植被在近红外“红边”区间的光谱特征信息，对监测植被长势与健康、水体Chl.a浓度、浮游植物分布等水生态信息灵敏有效[24]，运用红边波段的MCI与DoubleR两种模型可较好地检测内陆、沿海和海洋水域的表面水华和近地表植被信息[25-26]；而可见光区间内，绿光中心波长处Chl.a的光谱反射率，高于相邻两侧蓝光、红光波长的反射率，而形成典型的Chl.a反射峰，该反射峰强度(ρchl)与Chl.a浓度呈正相关[27-29]. 因此，本文基于Sentinel-2A影像获取Chl.a、TSM、MCI、DoubleR及ρchl5种指数，以此来探究漓江年内水质状况，具体计算公式为：

(6)

TSM=1.72IOPbp+3.1IOPbw

(7)

MCI=BVRE-BRed-0.53(BVRE-BRed)

(8)

(9)

(10)

式中，IOPapig表示浮游植物色素吸收特性，IOPbp表示海水中一般颗粒物的后向散射，IOPbw表示海水钙质白色粒子的后向散射，BVRE表示光谱中的红边波段，BRed表示光谱中的红波波段，BBlue表示光谱中的蓝波波段.

2.2 基于基本评价单元的漓江水文和水质信息定量分析

SDI反映岸线的不规则程度，刻画水面宽度差异，以此来研究漓江水体信息演变特征；WDr可定量表征基本评价单元内水域变化区域差异；土地利用动态度则适用于分析一定时间范围内某种土地利用类型的数量变化情况[30-32]，本文将土地利用动态度用于描述漓江水体动态变化分析，即计算水体信息变化动态度(K)，研究汛期、枯水期内漓江水文与水质信息动态变化程度. 因此，本文在278个基本评价单元内，采用SDI、WDr和K3个指标来定量分析漓江水体水文和水质的年内演变特征，具体计算公式为：

(11)

(12)

(13)

式中，S为水面周长；A为水域面积；Ua、Ub分别为基期、现期的数值；T为研究时段长，当T的时段设定为月时，K值就是该研究区某种水文信息的月变化率.

均方根误差(RMSE)[33]可反映观测值的可信度、精确度，本文利用均方根误差进行精度验证，计算公式为：

(14)

式中，xi表示观测值，yi表示真值，n表示观测次数.

3 研究结果

3.1 水体提取、叶绿素a和悬浮物质浓度的反演精度验证

为了探究NDWI、MNDWI、EWI、JQ、NDPI与后向散射系数S6种遥感指数提取漓江水体效果，本文选取漓江中、下游河段进行对比分析. 由图3可知，MNDWI指数能较好地提取出处于漓江中游的桂林市城区漓江水体，但在下游两岸山峰较多的区域，山体阴影和水体的区分不明显，MNDWI指数会把部分山体阴影错误归类为漓江水体，导致水域面积偏大. 相对于MNDWI指数而言，NDWI指数能较好划分水体和山体阴影，但在桂林城区段的水体提取效果不佳，临近建筑与水体难以通过阈值范围进行区分.EWI指数在山区的提取效果优于MNDWI指数，水体和山体阴影的区分比较明显，对城区的提取能力比NDWI指数强，能提取出较多的水体部分. 后向散射系数S完整提取整条漓江的能力比较强，可提取出河流水面宽度较小区域，但在下游会把部分山体阴影也划入漓江，整个范围内错误划分阈值的小图斑较多. 与其他指数相比，JQ与NDPI指数较好地区分出了桂林城区的水体，下游被错误地划分为水体的山体阴影面积均有所降低，细小河段以及河流连贯性也有所提高.

图3 基于Landsat 8、GF-1、Sentinel-1A与Sentinel-2A的6种指数最佳阈值提取漓江水体结果对比分析Fig.3 Comparison and analysis of the results of Lijiang River water extraction based on six index optimal thresholds of Landsat 8, GF-1, Sentinel-1A and Sentinel-2A

为验证主被动遥感加权指数JQ与NDPI两种水体指数的提取精度，选用汛期4月和8月与枯水期1月和12月水面提取结果进行对比分析，由图4中JQ指数与NDPI指数提取的水面对比结果可知，基于NDPI指数提取的漓江水体有较好的完整性、连贯性；表2为4个月水面宽度对比结果，其中JQ指数的水面宽度选取与实测水面宽度选取方式一致，NDPI指数的水面宽度则基于格网取均值；表3中NDPI指数提取水面宽度的RMSE在2.61～7.04 m之间，而JQ指数提取水面宽度的RMSE在6.05～9.79 m之间，对比两种数据提取误差范围，NDPI指数提取的数据整体误差明显小于JQ指数. 以上对比说明基于Sentinel-2A的NDPI指数提取水体信息有更高的精确度和可信度，下文将基于NDPI指数提取的水文信息分析漓江水体年内变化特征.

图4 JQ与NDPI指数水面提取结果对比Fig.4 Comparison of water surface extraction results of JQ and NDPI indexes

表2 JQ与NDPI指数提取水面宽度精度对比

表3 JQ与NDPI指数提取水面宽度的RMSE数值差异

利用36个Chl.a与TSM浓度站点监测数据，通过分析变化趋势、计算相关性来对比验证基于C2RCC算法的Chl.a、MCI、DoubleR及ρchl4种反演模型精度，并在每月上、中、下游反演数据中随机选取10个样点进行RMSE计算. 图5a～c中基于C2RCC算法的反演数据与实测Chl.a浓度变化趋势基本一致，使用SPSS统计分析软件进行相关性计算，在95%置信区间内R2为0.980，因每月河流径流量和水体流动速度不同，上、中、下游各部分Chl.a与TSM浓度不均，且河流两岸旅游开发程度与人类活动影响程度不同，所以通过随机选取的样点计算出的每月RMSE差异较大，结合表4可知1-12月Chl.a浓度的RMSE最小为0.18 mg/m3，最大为7.88 mg/m3；而图5d～f中MCI、DoubleR、ρchl与实测数据的归一化趋势可以看出，DoubleR的变化趋势吻合度较好，R2为0.906，MCI次之，R2为0.883，ρchl则不适于反映漓江Chl.a浓度变化，R2为-0.110；图5g～i中TSM浓度实测数据与反演数据变化趋势吻合度较好，R2为0.993，但上游12月与中、下游2月误差较大，1-12月中RMSE最小为0.17 g/m3，最大为12.55 g/m3，下文将基于C2RCC算法的反演数据分析漓江Chl.a与TSM浓度的年内动态变化规律.

图5 Chl.a、TSM浓度的遥感反演数据与实测数据的精度验证与趋势分析Fig.5 Accuracy verification and trend analysis of remote sensing inversion data and measured data of Chl.a and TSM concentrations

表4 Chl.a与TSM 浓度的RMSE差异

3.2 漓江流域水文与水质信息年内变化定量分析

利用从278个基本评价单元统计出的河流水面宽度、水域面积、Chl.a与TSM浓度，计算每月均值与标准差，并利用直方图(图6)进行趋势分析. 其中，均值表示一组数据的平均水平，反映数据集的集中趋势；标准差反映了数据集的离散程度；误差棒反应所测数据不确定度的大小.

由图6a和b可看出漓江1-12月水面宽度与水域面积的整体变化趋势基本相同，1-8月呈波动增长，8-12月呈波动型下降. 其中，3-8月上升趋势较明显，上升幅度也较大，3月开始进入雨季，降雨量增加，水面宽度、水域面积分别由3月的112.56 m、0.057 km2上升至8月的175.32 m、0.088 km2，最大值约为最小值的1.56倍；漓江降水主要集中在夏季，因此5-8月水面宽度、水域面积略高于其他月份. 图6c与d中Chl.a与TSM浓度变化趋势基本一致，2月富营养化程度最深，4-11月程度较低、水质较好，Chl.a与TSM浓度最大值分别为最小值的38与37倍；两者变化趋势与水面宽度、水域面积变化趋势有一定的负相关，即水面宽度、水域面积大的月份，水生植物、藻类等物质少，富营养化程度低，水质好. 为详细观察漓江水体Chl.a浓度的季节变化情况，选用离散度较低的1、4、7和10月进行分析，由图7可知每月Chl.a浓度最大值皆位于下游兴坪镇地区，该地区为阳朔县旅游重镇，旅游开发区较多，漓江水体富营养化最严重，4个月中Chl.a浓度最大值在7月，为39.27 mg/m3，而上、中游水体污染程度较高处均位于兴安县、灵川县与桂林市城镇居民区，其中1月上游水体污染区域较多，上游地区兴安县与灵川县居民区聚集、人口密度大，农田、厂矿分布集中，生活污水、工业废水等使得漓江水质变差.

图6 漓江1-12月水文与水质信息变化趋势Fig.6 Change trend of hydrological and water quality information of Lijiang River from January to December

图7 漓江水体1、4、7和10月各个河段Chl.a浓度对比分析Fig.7 Comparative analysis of Chl.a concentration in Lijiang River in January, April, July and October

为进一步探究漓江1-12月水体变化规律，本文将每月数据分为上、中和下游3部分，表5为水面宽度、水域面积、Chl.a与TSM浓度的全年上、中和下游均值与占比情况，图8描述了漓江逐月上、中、下游水体信息的动态变化趋势.

表5 漓江1-12月水文与水质信息上、中、下游均值与占比

根据表5可知，上、中、下游水面宽度、水域面积与Chl.a浓度变化较小，TSM浓度变化最大，且为中游>下游>上游. 结合图8可知，①上、中、下游水面宽度变化趋势相似，但存在差异. 1-3月水面宽度呈下降趋势，且2-3月下降显著，该期间气温回升但降水少，水面宽度明显偏低；3-8月雨量增多，水面宽度变大，且中、下游3-5月上升幅度最大，由114 m上升至178 m；8-12月水面宽度虽有波动但下降趋势明显，上、中、下游降幅分别为32、54和38 m，中游降幅最为显著；上、中游水面宽度在8月达到最大值，下游最大值出现在5月，上、中、下游最小值均出现在3月(图8a)，说明3月漓江流域径流量最小. ②水域面积与水面宽度总体变化趋势基本一致，说明水域面积与水面宽度有一定正相关性，即降水较多的月份，水面宽度增大，水域面积也增大(图8b). ③中、下游1-12月Chl.a浓度变化趋势相同，1-3月变化幅度较大，2月达全年最大值，说明该月份中、下游水体中水生植物最多，富营养化程度严重，而上游1-3月Chl.a浓度呈下降趋势；漓江水体4-10月Chl.a浓度均低于1.5 mg/m3，该时期降水较多，水体流动性大，水生植物的生长受到抑制，水质较好，而10-12月降水少，水生植物覆盖率随之上升，Chl.a浓度略有增加(图8c). ④TSM浓度变化趋势与Chl.a浓度基本相同，但上游1-3月变化趋势稍有差异，2月TSM浓度较低，该期间上游降水较多，水体流动性大，悬浮物覆盖率降低，该趋势与中、下游变化趋势相反(图8d). 由于1月数据离散度小，且基于基本评价单元提取的数据较完整，因此选取1月漓江流域278个基本评价单元的数据分析水文与水质信息的变化情况(图8e～h)，其中，水面宽度与水域面积变化趋势一致，上、中游水面宽度与水域面积略大于下游，两者下游变化趋势皆呈上升状态，但整条河流水面宽度与水域面积变化幅度较大，最大值分别为最小值的15.8和7.6倍；Chl.a与TSM浓度整体呈下降趋势，上游下降幅度最为显著，Chl.a浓度由第1个基本评价单元的11.99 mg/m3下降至第71个基本评价单元的2.63 mg/m3，TSM浓度由第1个基本评价单元的24.32 g/m3下降至第71个基本评价单元的0.59 g/m3，上游地区村庄较多，居民环保意识不足，致使河流污染严重；Chl.a与TSM浓度在中游第115个基本评价单元处有一个峰值，该区域的水域面积与水面宽度较小，表明1月桂林市该区域降水少、水生植物与悬浮物质较多.

图8 漓江水文与水质信息变化趋势Fig.8 Change trend of hydrological and water quality information of Lijiang River

3.3 漓江流域水文与水质信息变化差异的定量分析

水面变化区域差异值可以比较直接地观察河流水面变化的情况，WDr值越大说明变化程度越大. 以全年数据均值作为Ua，每月均值作为Ub，计算河流信息变化差异值，并以折线图来刻画每个月份的变化情况.

由图9可以看出，漓江水面宽度、水域面积、Chl.a和TSM浓度变化与图8中变化趋势呈正相关. 图9a和b中水面宽度WDr与水域面积WDr变化趋势相同，4条差异变化折线在8-10月都是正值，即为正向变化，3、4和12月都为负向变化，3月变化差异最大，水面宽度与水域面积流域WDr分别为-0.24与-0.23. 漓江上游段在3-7、12月为负值，即相对全年来说，该期间降水较少，水面宽度和水域面积处于减小状态，而1-2、8-11月降水较多. 其中，3月干旱严重，负向变化差异最大，水面宽度差异为-0.28，水域面积差异为-0.30. 漓江水面宽度与水域面积的流域变化差异与中、下游部分的变化趋势大体一致，且中、下游5-10月雨量大，水面宽度、水域面积较大. 中、下游正向变化差异最大的月份为5和8月，负向变化最大的仍为3月. 漓江3-12月各河段Chl.a浓度变化趋势相同，4-11月WDr皆为负值，说明该时期Chl.a低于全年平均浓度，而1-3月上游Chl.a浓度变化趋势与中、下游部分稍有差异，该时期上游WDr值虽呈下降趋势，但Chl.a均高于全年浓度(图9c)；漓江中、下游Chl.a浓度变化与整条流域变化趋势相同，2月正向变化差异最大，总流域WDr为4.16，即与全年相比2月Chl.a浓度增幅最大，同时说明该月份水生植物覆盖率最高. 漓江TSM与Chl.a浓度变化趋势大体一致，两者变化趋势与水域面积、水面宽度的变化趋势呈一定的负相关；上游1-3月TSM浓度与中、下游变化趋势相反，2月TSM浓度增幅最小(图9d). 1月漓江水面宽度WDr与水域面积WDr整体变化趋势较稳定，但仍有部分区域变化幅度较大；水面宽度与水域面积下游区域WDr负值较多，两者负值基本评价单元数量分别占总负值数量的45.14%与46.56%，即与全年相比，1月下游地区降水较少，水面宽度、水域面积多处于减小状态(图9e，f). 而漓江水体Chl.a的WDr与TSM的WDr皆呈下降趋势，其中上游Chl.a的WDr均为正值，TSM的WDr仅有第26个基本评价单元为负值，上游兴安县与灵川县内漓江水质较差；中游42%格网Chl.a的WDr为负值，85%格网TSM的WDr为负值，其中第113个格网正向变化最大，WDr达到了1.96；下游81.48%的格网Chl.a的WDr为负值，TSM的WDr均为负值(图9g，h)，说明1月漓江上游水面水生植物与藻类等悬浮物质覆盖率最大，富营养化程度最深，下游水体水面宽度与水域面积虽呈减少趋势，但水质较好.

图9 漓江水文与水质信息变化区域差异Fig.9 Regional differences of hydrological and water quality information changes of Lijiang River

3.4 漓江流域河岸线发育系数定量分析

河岸线发育系数可表征岸线的不规则程度，也可表征水面变化程度. 漓江河流SDI与水面周长变化呈正相关，而与水域面积变化呈负相关(图10). 由图10a和b可知，1-12月水面周长与SDI皆有下降趋势，但岸线发育系数变化幅度较大，可分成3个阶段：小幅下降(1-3月)、急剧下降(3-5月)和稳定(5-12月). 1-3月雨量较少，流域水量低，水域面积小，河岸线较为曲折，河流水面周长较大，但仍有下降趋势. 3-5月降水增多，河流面积迅速增加，该期间河岸线变的平缓，河流水面周长降低. 5-12月SDI值在1.6～1.7之间波动，岸线发育率较低，该阶段河流水面周长虽有下降趋势，但由于雨期的来临，水量增多，水域面积明显大于1-4月，岸线不规则程度也较低. 图10c表征了漓江1月278个基本评价单元的SDI变化情况，整体来看该月份岸线发育率呈下降趋势，但中游起伏较大，其中，第106个基本评价单元的SDI值达到5.89，桂林市内该水域面积数值较小，水面周长达到最大值，说明此基本评价单元中的水量少，岸线较为曲折.

图10 漓江河流水面周长和岸线发育系数的变化趋势Fig.10 Variation trend of river circumference and coastline development coefficient of Lijiang River

表6 基于汛期和枯水期的水文与 水质信息动态度变化

3.5 漓江流域水文与水质信息汛期、枯水期变化定量分析

为探究漓江汛期和枯水期水体信息动态变化规律，本文以枯水期均值作为动态度公式中的Ua，汛期均值作为Ub，来计算汛期相对于枯水期的变化情况. 由表6可知，水面宽度与水域面积的总体流域动态度都是正值，Chl.a与TSM的总河流动态度都是负值，即汛期时降雨多，漓江水面宽度、水域面积都大于枯水期，而水体中的Chl.a与TSM浓度低于枯水期，其中Chl.a与TSM浓度的变化较大，动态度均小于-7%，下游水面宽度与水域面积动态度分别为2.24%与2.41%，说明水中植物与TSM浓度变化最大；上游水面宽度与水域面积为负值，即上游汛期水面宽度与水域面积小于枯水期，说明上游汛期降水少于枯水期；而中、下游降水多于枯水期，且下游降水最多、动态变化最大；汛期时漓江水质好，富营养化程度低，整条河流Chl.a与TSM浓度均小于枯水期，其中下游变化最大，动态度为-13.09%.

本文基于278个基本评价单元进一步分析各河段漓江汛期和枯水期水体信息动态变化，具体见图11. 由图11a可知，汛期大部分地区的水面宽度在100～250 m范围内波动，但上游部分地区变化趋势较大，其中最大值为灵川县第31个格网，水面宽度达到了338.31 m，最小值为第5个格网位于兴安县境内，数值为40.97 m. 枯水期上游水面宽度波动较大，中游水面宽度有下降趋势，下游略呈上升状态(图11b). 汛期与枯水期水域面积与水面宽度变化趋势基本一致(图11c，d)，但中游第105个基本评价单元水域面积较大，汛期与枯水期水域面积分别为0.14和0.17 km2，该地区全年降水多、水量多、水面宽. 汛期与枯水期Chl.a浓度整体呈“双谷型”变化(图11e，f)，中游第100～110个基本评价单元处浓度较高，该地区村庄聚集，较多生活污水的排放导致河流富营养化程度加深. 除第26个基本评价单元TSM浓度为16.85 g/m3，达整条河流的最大值外，汛期TSM浓度均在0～3 g/m3之间(图11g)，该水域的水面宽度、水域面积与Chl.a浓度也较大，说明汛期该地区水量多，但水质差. 枯水期TSM与Chl.a浓度呈正相关(图11h)，上游兴安县第1～10个、中游桂林市区第10～110个基本评价单元所在水域水面富营养化程度高，其中第107个基本评价单元处TSM浓度为15.03 g/m3，达整条河流的最大值，该水域的水面宽度、水域面积与Chl.a浓度也较大，说明该处水体污染程度最严重.

图11 漓江汛期和枯水期水文与水质信息变化趋势Fig.11 Change trend of hydrological and water quality information in flood and dry seasons of Lijiang River

4 结论与讨论

本文以2016-2020年Landsat 8 OLI、GF-1和Sentinel-2A多光谱数据和2018年逐月Sentinel-1A IW模式的SAR影像为主要数据源，利用构建的多种主被动遥感水体指数提取了漓江逐月的水域信息，选用4种模型反演了漓江逐月的Chl.a与TSM浓度，并基于基本评价单元采用多种水体指标定量分析漓江水文、水质的年内、枯水期和汛期上、中和下游的时空演化特征，结果表明：(1)主被动遥感加权指数JQ与NDPI指数的提取效果优于NDWI、MNDWI、EWI与后向散射系数，而NDPI指数与JQ指数相比提取精度更高、可信度更强，水面宽度RMSE在2.61～7.04 m之间. (2)基于C2RCC算法反演的数据可较好地揭示漓江Chl.a与TSM浓度的年内动态变化规律，Chl.a浓度的RMSE处于0.18～7.88 mg/m3之间，R2为0.980，TSM浓度的RMSE为0.17～12.55 g/m3，R2为0.993，而DoubleR可较好地反映Chl.a浓度的变化趋势，R2为0.906，MCI的R2为0.883，ρchl的R2为-0.110，则不适用于反映水体Chl.a浓度变化. (3)1月水面宽度与水域面积变化较稳定，大部分基本评价单元所在水域的水面宽度和水域面积分别处于60～240 m与0.03～0.12 km2区间内，而Chl.a与TSM浓度明显呈下降趋势，其中上游降幅最大，由第1个基本评价单元至第71个基本评价单元Chl.a与TSM浓度降幅分别为9.36 mg/m3与23.73 g/m3，中、下游变化则较稳定，实测数据则依靠站点监测，所得数据较分散，无法进行连续性分析. (4)漓江1-12月平均水面宽度在100～175 m范围内，水域面积处于0.05～0.09 km2之间，其中 5-10月降水较多，水面宽度在150 m以上，水域面积在0.07 km2以上，水面宽度和水域面积最大值出现在8月；而5--10月水体富营养化程度低，水质较好，Chl.a浓度低于0.6 mg/m3，TSM浓度低于1.3 g/m3，漓江2月水质最差，Chl.a浓度为9.82 mg/m3，TSM浓度为15.32 g/m3，水体富营养化程度较高地区主要集中于上、中游的兴安县、灵川县等城镇居民区，以及下游旅游开发较多的兴坪镇. ②漓江水体水面宽度WDr与水域面积WDr在8-10月都大于0，即为正向变化，3、4、12月全为负向变化，其中3月负向变化最大，水面宽度与水域面积的WDr分别为-0.24与-0.23，该月份降水最少；漓江流域4-11月Chl.a的WDr均小于0，该时期Chl.a浓度低于全年平均值，TSM浓度的WDr值也在0附近波动，该时期漓江水质较好，富营养化程度低. 漓江1-12月水面周长与SDI呈下降趋势，4-5月降幅较大，其中上游降幅为0.57，下游降幅为0.27，5-12月SDI则基本稳定在1.6～1.8之间，该期间水面变化程度较小. 漓江水体水面宽度、水域面积与Chl.a浓度在上、中、下游总体变化较小，TSM浓度变化最大，且为中游>下游>上游；汛期时降水多、水体流动性强，大部分地区平均水面宽度在100～250 m范围内，但上游地区波动较大，TSM浓度则基本处于0～3 g/m3之间，水体富营养化程低，水质较好.

本文利用Landsat 8 OLI、GF-1、Sentinel-1A与Sentinel-2A 4种卫星影像提取漓江水面数据，并选用C2RCC算法、MCI、DoubleR与ρchl进行漓江水质反演，研究结果与漓江实测数据较为吻合，但还存在一些不足，将会在以后工作研究中逐一完善：1)因部分月份中影像受云雾和阴雨天气影响，部分星地数据存在不同步现象. 2)在C2RCC处理过程中，所有参数均依据美国宇航局生物光学海洋算法数据集的训练结果进行设置，但为进一步提高实验结果的准确性，应尽可能根据全球各地不同近岸水体特点，通过当地水质实测数据来优化调整为更有区域适应性的水质参数. 3)主被动遥感加权指数JQ中各指数的权重是通过试错法反复试验所设定的，存在主观赋值所带来的结果不确定性等问题，同时还应当考虑不同河段地物类型的特点，进行分段差别化处理.

5 附录

附表Ⅰ见电子版(DOI: 10.18307/2021.0306).