高分一号WFV影像的深圳市水库CODMn浓度时空分布特征

李俊，张文志，邓孺孺，鲁志文，梁业恒，沈雪娇，熊龙海，刘永明

1.中山大学 地理科学与规划学院,广州 510275;2.广东省水环境遥感监测工程技术研究中心,广州 510275;3.广东省城市化与地理环境空间模拟重点实验室,广州 510275;4.广东省水文局惠州水文分局,惠州 516003

1 引 言

随着经济快速发展、城市化进程不断加快和人口数量的加剧，对水资源的依赖逐渐加重，增加了对地表水的需求和地表水污染的风险。城市及其附近的地表水以有机污染为主，了解地表水的有机污染状况和变化趋势可以有效帮助保护紧缺的水资源，为水环境治理提供参考。

传统的点采样监测方式虽然可以给出准确的测量结果，但是既耗时又花费高，并且无法提供大范围面状的水质评估和实时的水质状况。遥感技术凭借探测范围广、成本低、周期性重复覆盖等优势，在水体有机污染监测中广泛使用（Gholizadeh 等，2016；Kutser，2012；Prangsma 和Roozekrans，1989）。

国内外学者都对水质遥感中的有机污染情况进行研究。国外，Gitelson 等（1993）建立了一套光谱反射率与水质参数的经验关系，同时反演内陆水体的叶绿素、悬浮物、CDOM，并在巴拉顿湖得到应用；Vignolo 等（2006）构建水质参数WQI模型追踪淡水系统中的有机污染；Alparslan 等（2009）利用Landsat 5 TM 数据，建立了土耳其的Kucukcekmece 湖卫星反射率与COD、叶绿素a 等水质参数的多元回归关系；Fichot 等（2016）使用高分辨率成像光谱监测旧金山湾三角洲河口重要水质指标，包括浊度、悬浮颗粒、CDOM、叶绿素等。国内，傅江等（1993）利用野外实测光谱及同步的实测数据建立苏南大运河COD 含量的线性回归模型；Wang 等（2004a）结合TM 影像和现场测量数据反演深圳市水库生物需氧量和化学耗氧量；Chen 等（2007）用SeaWiFS 离水反射率数据反演珠江口水体综合污染指数，其结果与WQL 图具有相似分布；Wang 等（2008）基于神经网络技术利用TM 数据反演内陆水体的多项水质因子；黄妙芬等（2011）通过反演水体CDOM 光学量，基于相关性分析估算辽东湾水体COD 含量；Li 等（2016）使用融合高空间分辨率影像和高光谱影像的经验算法来监测水体溶解氧、高锰酸盐指数（CODMn）和氨浓度。目前的有机污染遥感研究主要基于经验模型（陈玲侠，2017；郭丽峰等，2010；李佳琦等，2019；张雪等，2019），依赖覆盖各种环境条件的大量实测数据，受到时间和地区的限制；有机污染程度评价主要利用污染指数、化学需氧量、CDOM 光学量（Li 等，2017；Shirke 等，2016；邓孺孺 等，2002）等水质参数不明确的综合指标表示。

基于此，本文以深圳市3 大水库为研究对象，从具有更强稳健性和普适性的辐射传输模型（Stumpf 等，2016）出发，根据光在水体内部的辐射传输过程，通过测量获得研究区水体固有光学参数，建立反射率与水质参数CODMn之间的物理模型。选用GF-1 WFV 数据，反演3 大水库的CODMn浓度，分析CODMn浓度时空变化特征，提供水质保护应对策略分析，为深圳市饮用水源地水环境保护提供参考。

2 研究区与数据

2.1 研究区概况

深圳市位于东江流域下游，75%的饮用水源来自东江，主要的调蓄供水水库有铁岗水库、西丽水库、石岩水库等。3 大水库位于东部引水工程末端（图1（a）），库容量分别为深圳市第1、第3和第4，水库蓄水量大，可以很好代表深圳市水库水质状况。石岩水库有6条入库支流，与石岩街道相邻，南部建有铁岗水库—石岩水库供水渠与铁岗水库相连。铁岗水库有4条入库支流，东部有西铁连通隧道与西丽水库相连。西丽水库主要调蓄深圳东部水源至铁岗水库。3 大水库紧靠深圳市建成区，库区周边分布一些小型加工厂和住宅，周围用地大多为建筑用地、果园和菜园，外围用地基本是经济林、农田和水塘。3大水库存在潜在的污染源威胁库区水质安全（图1（b）），水库基本概况见表1。

表1 石岩水库、铁岗水库、西丽水库概况Table 1 Features of Shiyan reservoir，Tiegang reservoir and Xili reservoir

图1 深圳市3大水库位置图Fig.1 Location of Shenzhen’s three major reservoirs

2.2 现场水样数据

CODMn现场样点是采集表层（0—50 cm）水样，冷藏保存，同时利用GPS 记录每个采样点坐标。在采样完毕后迅速送往实验室，采用GB 11901-89测定水样中CODMn浓度。其测定原理为采用硫酸—高锰酸钾消解体系，水样经沸水浴加热消解后，加入过量草酸钠，再用高锰酸钾标准溶液回滴，计算出高锰酸盐指数。

2.3 遥感数据

本研究选择太阳同步回归轨道的GF-1 卫星16 m 空间分辨率的多光谱CCD 遥感数据，幅宽800 km，重访时间2 d。传感器为WFV3，有4个波段，波段设置分别为蓝光波段（450—520 nm）、绿光波段（520—590 nm）、红光波段（630—690 nm），近红外波段（770—890 nm），中心波长分别为514 nm、546 nm、656 nm 和822 nm。选择成像时间从2018-03-11—2019-05-12 的遥感影像进行处理。

3 CODMn遥感反演原理

3.1 水体光谱特征

考虑到3大水库与广东省内陆水体水质因子物质来源类似，在广东省内陆水域选择典型单一类型主导的水体样品进行测量，可认为其隔离了其他物质干扰，典型水体样品水质参数见表2，其平均离水反射率的曲线见图2。

表2 野外典型水体样品水质参数Table 2 Water quality parameters of typical water samples

图2 典型水体样品离水反射率Fig.2 Water reflectance of typical water samples

耗氧性有机物对水体光谱的影响主要表现为吸收性，与叶绿素和悬浮泥沙对水体反射率的增强作用相反。耗氧性有机物的吸收包括溶解性有机物的吸收和有机碎屑颗粒物的吸收，对可见光的吸收作用很强。随着水体有机污染程度增加，水体可见光范围内的离水反射率变低，水体颜色加深，达到蓝黑、甚至黑色。可见光光谱内，对绿光、红光波段影响最为明显，降幅最大（王云鹏等，2001）。在420—520 nm波段范围内，CODMn浓度与水体反射率呈负相关（傅江等，1993）；在550—580 nm 波段范围内，受到有机污染水体反射率出现峰值，而波峰较正常水体平缓；在700 nm的反射率变大，在污染程度较大时，不具备明显反射峰（温爽等，2018；赵起超等，2019）。耗氧性有机物的存在明显影响到水体光谱特征。

3.2 CODMn指数

CODMn指数是指在一定条件下，用高锰酸钾氧化水样中耗氧性有机物和还原性无机物所消耗的氧量。在水体交换能力强的地表水中，由于表层水与大气氧直接接触，致使大多数还原性无机物难以留存，CODMn指数常用来反映的是水体被耗氧性有机物污染的程度（陈静生等，2003）。

在中国，CODMn、溶解氧、生物需氧量等水质参数被用来反映水体有机污染程度。由于水体中溶解氧易受到温度的显著影响，生物需氧量受细菌种类和其他因素的影响，CODMn更能反映水体有机污染程度（Xia 等，2005）。对于有机污染物含量较低的水体，适合用高锰酸盐指数（CODMn）法进行评价（Vaillant 等，2002）。CODMn作为一个综合指标，评价内容不仅包含可溶性有机物，还包含不溶于水的颗粒有机物和胶状腐殖质。当CODMn值超过4 mg/L 时，表明水体已经受到有机物的污染（林碧琴和谢淑琪，1988）。

3.3 CODMn遥感反演模型

分析太阳光在水体的辐射传输过程，到达水面的电磁波一部分发生镜面反射进入大气，一部分传入水体。传入水体的电磁波与水体及水体中的组分发生吸收散射作用，在光学浅水区，电磁波还经由水底部分反射再进入水体。离水光谱特征受水体表面反射与折射、整层水体吸收散射、底质吸收反射作用等影响。卫星传感器接收到的离水辐亮度包含丰富的水体杂质的类型和浓度、底质种类、水深等信息，为水质遥感反演提供理论基础。

刚好在水面以下的上行辐亮度（Lu-）包括水体上行散射的贡献（Lwc）和水体底质部分的贡献（Lbot）（Philpot，1989）：

根据比尔定律，辐射透过总吸收系数为a，散射系数为b的均匀水体，在水深为h的下行辐照度：

式中，Ed-为刚好在水面以下的下行辐照度。水体散射部分以薄层水分析，如果只考虑一次散射，在水深为h处，厚度为dh的薄层水的散射作用：

式中，为水体散射相函数，对式（3）积分，令水体总消光系数c=a+b，则深度为h的整层水体上行散射部分的辐亮度贡献：

假设底质平坦且符合朗伯体辐射规律，底质反射率为ρ，则水体底质部分的辐亮度贡献：

以反射率表示辐射传输过程的关系：

水体光学特征主要由水体本身的性质决定，同时又受到各种水质因子的影响。在水环境遥感研究中，考虑的水质影响因素以悬浮泥沙、叶绿素、CDOM 为主。CDOM 仅占水体溶解有机物总量的50—70%（AiKen 等，1985），并且广东省内陆水体特征与大洋水体不同，有机污染物除了溶解有机物还含较多的碎屑、胶状腐殖质、不溶于水的颗粒有机物。颗粒有机物与可溶有机物的吸收光谱较为相似，无法从光谱特征有效的区分开来（秦雁，2015）。CODMn作为综合水质参数，更适宜评价研究区有机污染情况。通过对典型重CODMn水体光学参数准确测量，用CODMn表征水体耗氧性有机物，将CODMn加入到影响水色的要素中，构建CODMn遥感反演模型。

对于水体综合光学参数，式中下标w、s、chl、p分别表示纯水、悬浮泥沙、叶绿素、耗氧性有机物，D为各组分物质的对应浓度。

基于反射率的关系式（6），水体各组分光学参数可以通过测量获得，模型中的未知数仅为研究区水色3种组分浓度Ds、Dchl、Dp、水深h、底质反射率ρ。对于未知数多于波段数的病态方程组，需分析研究区水体光学特性，选择敏感波段建立方程组，即可求解研究区水体的CODMn浓度。

4 深圳水库CODMn遥感反演

4.1 光学参数测量

CODMn遥感反演模型的水体固定光学参数需要测量4类物质的光学参数。纯水吸收系数和散射系数参考（邓孺孺等，2012a）方法进行实验。叶绿素、悬浮泥沙、耗氧性有机物吸收系数和散射系数测量参考（邓孺孺等，2015）方法，使用实验室自主研制的消光系数测量装置分别测量表1典型水体样品的消光系数，根据光学深水中近似认为h→∞，计算得到散射系数和吸收系数光谱。其中，耗氧性有机物吸收系数、散射系数光谱见图3。根据GF1-WFV 波段设置，计算得到对应波段的吸收和散射系数。

图3 耗氧性有机物吸收、散射系数Fig.3 Absorption and scattering coefficient of oxygen-consuming organic

4.2 数据预处理

GF-1 WFV3 数据预处理主要包括几何纠正、辐射定标、大气纠正、水陆分离。

（1）几何纠正以高空间分辨率的Google Earth影像为基准，选择的GCP 点均匀分布在影像内部。通过控制点的选择，确保校正误差小于0.5 个像元。

（2）高分一号辐射定标根据中国资源卫星中心官网（http://www.cresda.cn［2019-11-14］）提供的2018年国产陆地观测卫星绝对辐射定标系数，结合头文件提供的卫星图像捕获参数，将DN值转换为行星反射率。

式中，ρTOA是行星反射率，d是卫星过境时的日地距离和平均日地距离的比值，F0是大气层外太阳辐照度。

（3）大气校正是准确获取水色参数的前提条件（马荣华等，2009），本文针对高分一号波段设置，充分考虑天空光的下行散射影响。采用暗像元法，选择图幅内水质级别长期为Ⅰ级的开阔的清、深水体作为暗像元进行大气校正。

假设天空辐照度各向同性，地表是朗伯面反射，忽略大气分子和气溶胶的相互作用，天空光部分仅考虑瑞利散射作用。建立行星反射率与真实地物反射率（ρs）与的函数关系，形式如下（齐志新等，2007）：

式中，大气总散射系数ω由暗像元算出，忽略大气的吸收作用，观测方向的大气透过率tu= 1 -ω，根据大气光学厚度定义求出入射方向的大气透过率td。为提高大气校正精度，考虑大气的非均质性，选择多个清、深水体像元，根据图幅像元到各暗像元距离的权重，计算出各像元对应的大气校正参数，根据式（10）将行星反射率转换为地物反射率。在大气校正的同时，参考邓孺孺等（2012b）方法去除水面镜面反射的信息。从影像上分别选择水体和植被两种典型地物样本，对比大气校正前后纯像元反射率和实测地物反射率，对比结果如图4 所示。由图4 可以看出，大气校正消除大气散射光对蓝、绿、红波段反射率的增强，减弱大气散射光对近红外波段反射率的衰减。经过大气校正后，两种典型地物的反射率与实测地物反射率很接近。

图4 大气校正前后水体、植被反射率同实测值的比较Fig.4 The reflectivity before and after the atmospheric correction，and measured values of vegetation and water

（4）水陆分离，提取水体部分，能够突出水体部分的信息，有效提高运算速度。基于光谱特征的水体提取中，单波段或多波段阈值法简单易用，由于GF-1 卫星WFV 数据只有4 个波段，综合考虑，决定选用水体指数（NDWI）方法提取水体。NDWI方法是目前广泛应用的一种水体监测方法，利用水体在绿光波段的高反射率和近红外波段的低反射率特征，增强水体信息；利用植被、土壤在近红外波段的高反射率特征，抑制植被、土壤信息。NDWI 同时还可以提供水体浊度信息（Mcfeeters，1996），在二类水体遥感监测中，能有效增加水体提取精度。本文利用近红外波段和绿波段中像元亮度值进行计算，其定义：

结合实际中存在浑浊度较高的水体，依据水体近红外波段反射率大于绿光波段的特征，加入NIR设置水陆分割阈值。图5 显示，深圳市3 大水库细节部分能成功提取，保证了水陆分离精度。

图5 2019年1月25日影像水陆分离结果图Fig.5 Land and water separation results of image on January 25，2019

4.3 特征波段选择

从水体弱信息中提取CODMn的关键是利用敏感波段区分耗氧性有机物与水体其他组分。针对GF1 WFV 的波段设置，近红外波段对应耗氧性有机物散射太弱，水体光谱主要受到悬浮泥沙、水体自身吸收等的影响，难以提取CODMn微弱的信息，故排除。绿光波段对应耗氧性有机物有反射峰值，对应叶绿素和悬浮泥沙均表现为弱吸收、弱散射，将其列入特征波段；蓝光和红光波段可以提供足够信息提取水体有机物含量（Bowers 等，2004），故将蓝光和红光波段加入特征波段。综上，最终确定CODMn特征波段是蓝、绿、红3 个波段。

4.4 CODMn反演

数据预处理后，得到遥感反射率影像。通过下式（Lee 等，2002）转换为水面以下的反射率。

考虑深圳地区水库水深较大，可以认定为光学深水，近似认为h→∞，忽略底质影响。根据研究区光学深水的特性，简化CODMn遥感反演模型得到：

水体散射相函数按照瑞利散射的散射相函数进行计算。将实验测量的纯水、耗氧性有机物、叶绿素、悬浮泥沙光学参数数值参考式（7）和式（8），结合蓝、绿、红共3个波段代入到式（13）中建立3 个方程组，联立求解水体综合光学参数，来得到整幅影像的CODMn浓度Dp。影像CODMn浓度提取图进行密度分割将浓度分级，得到图6。

5 结果与分析

在一定区域内物质来源稳定的情况下，耗氧性有机物的光学参数也是稳定的，本文基于广东省内陆水域反演结果来评价CODMn反演模型的精度和可靠度。由于水质污染情况稳定，日期间隔较短，本次研究验证数据来自广东省水文局惠州水文分局提供的准同步的监测站数据和取水样分析数 据。对 比2018年3月11日、2019年1月25日CODMn浓度模型反演结果与实测真实值，两次模型验证点分别为15 个和21 个，验证点分布位置如图6所示，验证点实测值与模型反演值见表3。

图6 CODMn浓度反演结果与验证点分布Fig.6 CODMn concentration inversion results and distribution of verification points

表3 实测值与模型值统计对比Table 3 Statistical comparison of measured values and model values

采用统计指标决定系数R2、相对误差RE 和均方根误差RMSE作为检验模型精度的标准。

实测值与模型反演结果线性关系明显（图7），决定系数R2分别为0.832 和0.815，RMSE 均方根误差分别为0.464 mg/L 和2.220 mg/L。2019年的验证中，17 号点的误差最大，远高于实测值结果。经实地考察，17 号点汇入支流较多，且支流水质较差，有水质不稳定的情况，如图8所示。从模型精度评价来看，模型应用在广东省内陆水体的多个时相，CODMn反演结果的空间规律与实际规律一致，且反演精度较优，说明基于广东省内陆水体的光学特征，从辐射传输机理出发建立的CODMn遥感反演模型能很好地应用在具有相似物质来源的水域中。验证结果表明本文建立的遥感反演模型方法可靠，能应用在深圳市3 大水库的CODMn浓度反演，结果有效。

图7 模型反演值与实测值对比Fig.7 Comparison of model inversion and measured values

图8 17、18号遥感反演结果与实地考察照片对比Fig.8 Comparison of remote sensing inversion results with field investigation photos of No.17 and No.18

5.1 空间分布特征

图9显示，深圳3大水库CODMn整体浓度不高，受到轻度有机污染，仅在库角处出现浓度高值，高值区浓度一般为14.1—17.6 mg/L，峰值可以达到24.0 mg/L。3大水库整体。石岩水库的库东南与石岩河相接靠近居民区的位置出现大范围有机污染，由于库坝隔离，CODMn浓度出现明显线性隔断（图9（a）、（e））。铁岗水库多处库角都显示CODMn浓度偏高，主要是库西南角、库北、库西北，虽然库北与石岩水库相通，但由石岩水库连接处水质推断，受到有机污染水体来源并非石岩水库。西丽水库主要是库东和库北处CODMn浓度高，分别位于与东江和白芒河相连接处。深圳3 大水库中，污染面积较大的是石岩水库，铁岗水库和西丽水库水质污染程度情况相对较低。

图9 2018年3月—2019年5月CODMn浓度反演结果Fig.9 CODMn concentration distribution from Mar 2018 to Jan 2019

表4 模型反演值与实测值的误差Table 4 Errors between model inversion and measured values

深圳3 大水库CODMn浓度空间分布特征是总体有机污染程度较轻，但在局部库角出现CODMn浓度高值，主要是分布在与河流连接处。结合水库周边土地利用类型，受污染处多接近生活区，临近处有建筑物，推测主要受人类活动影响。

5.2 时间分布特征

选 择2018年3月11日—2019年5月12日 共计6 景影像，包括4 个季节。根据深圳市气象特征划分，4—9月为丰水期，10—11月为平水期，12—3月为枯水期（宋云龙等，2017）。3大水库CODMn浓度峰值和平均值的变化如表5所示，最小的浓度峰值出现在2018年7月。在丰水期，3 大水库CODMn平均浓度和库角CODMn浓度都明显低于枯水期和平水期，且有机污染范围减小。结合气象监测结果，由于前一天的降雨，库区水域面积增大，降雨对污染物的稀释作用明显。四季中，3 大水库CODMn浓度平均值在11月最高，在7月和1月最低，秋冬季节CODMn浓度峰值明显高于夏季。2018年3月至11月，石岩—铁岗水库耗氧性有机物有向南转移的趋势，但库区CODMn浓度平均值呈降低的趋势。从年度变化来看，2018年—2019年3 大水库有机污染总体减弱，水质得到改善，至2019年5月CODMn平均浓度达到2.27 mg/L，浓度分布较均匀，无明显浓度差异，库角处几乎无浓度峰值出现。

表5 深圳市3大水库CODMn浓度峰值和平均值Table 5 Maximum and average CODMn concentrations of the three major reservoirs in Shenzhen

5.3 水质保护应对策略分析

石岩、铁岗、西丽水库水体均位于饮用水水源保护区内，水质保护效果较好，据2018年度深圳市环境状况公报，西丽水库、铁岗水库水质达到国家地表水Ⅱ类标准，石岩水库水质达到国家地表水Ⅲ类标准。但深圳的饮用水水源地仍存在用地规划靠近水源保护区，水源保护区内留有大量建筑物的现象，如图10（a）所示。

图10 石岩水库保护区违建拆除前后对比Fig.10 The demolition of illegal buildings in Shiyan Reservoir Protection Area

2018年，深圳市政府进行年度治水提质大会战，解决全市集中式饮用水水源地89 项问题，拆除水库保护区违建、关停保护区工业企业、整治面源污染。其中石岩水库保护区部分违建拆除情况如图10 所示。经过一年的专项整治，2018年3月—2019年5月深圳市3 大水库水质出现好转，如图9所示。

3 大水库邻近陆地土地利用类型多样，有工业用地、居民用地、农业用地等，受人类活动影响不可避免。同时，蓄积在库内的内源污染会在无外源污染输入及上层水质较好的情况下爆发水质污染（许士国和汪天祥，2015）。水库底泥中有机质、氮和磷含量较高（聂祥，2009；袁文权，2004），可以在短时间内消耗上覆水体的溶解氧，造成水体CODMn浓度急剧变化。针对内源污染状况，可以采取疏浚底泥、修建人工湿地，加强水库水体生态环境的自我修复能力。结合遥感影像，水库周边面源污染源供给和人类活动对深圳市3大水库CODMn有着决定性贡献，3 大水库水质保护的核心是控制外源污染体，特别是丰水期，避免污染源的输入。

6 结 论

CODMn作为《国家地表水监测及评价方案》中9 大基本监测指标之一，直接参与水质评价。CODMn浓度的准确和长期监测对水环境保护和治理具有重要意义。针对目前CODMn反演局限在依赖大量实测数据的经验模型，本文基于辐射传输原理，考虑广东省内陆水体特征，将影响研究区水色反演的主要水质组分：悬浮泥沙、叶绿素、耗氧性有机物考虑到模型中，测量耗氧性有机物光学参数，建立地物反射率与CODMn浓度之间模型。通过对深圳市3 大水库CODMn浓度进行定量反演，经过地面测量数据验证，决定系数R2达到0.832，均方根误差较低，模型方法可靠。本文模型方法与线性回归法、波段比值法等经验模型相比，优点在于模型中参数意义明确，反演机理清晰，计算相对简单，可以通过改变不同地区水体光学固有参数来适应不同的地区的水体CODMn浓度反演。

对深圳3 大水库CODMn浓度时空分布进行研究，主要得出以下结论：深圳3大水库的水质整体受到轻度有机污染，CODMn高值区多位于库角处。水库间连接处未出现污染转移扩散情况。随着时间推进，3 大水库水质情况得到改善。结合2019年5月份水质结果，验证深圳市政府2018年度水质整治专项措施有效性。经清理水库周边违建、转移居民等切断潜在外源污染体后，3 大水库水质得到改善。水库水质保护的核心是控制外源污染，避免丰水期污染源的输入。

本文在进行CODMn浓度反演时，仅依赖一组水体固有光学参数，而水体固有光学特性存在着季节差异，可能会产生一定的误差，后续将充分研究水体固有光学参数的变化模式，以提高模型的稳健性。另外，本文模型应用在内陆浅水时，引入水底反射的影响，造成模型未知变量过多，需要更多波段的遥感数据，有待开展更多谱段的卫星数据如高分六号在浅水CODMn反演方面的研究。