基于主波长水色参数的水库采砂遥感监测

秦 雁，吉红香，黄本胜，陈亮雄，周 宇，杨静学

(1.广东省水利水电科学研究院，广州 510635； 2. 河口水利技术国家地方联合工程实验室，广州 510635； 3.广东省粤港澳大湾区水安全保障工程技术研究中心，广州 510635)

1 研究背景

采砂是一项剧烈干预河道自然演变的经济活动[1]。适当有序的采砂不仅有利于河道行洪及航运通航，且能提供宝贵的砂石资源[1-2]。近年来，我国城市建设与基础设施建设等经济活动对砂石需求量持续增长，高利益驱使下各种非法采砂活动屡禁不止。盲目、超量和无序采砂对防洪、涉水工程、供水、航运、渔业、水生态安全以及社会稳定产生恶劣的负面影响[1-5]。然而，目前采砂监管工作主要是以执法船只或车辆巡逻的方式进行的，面临着监管区域广、专职人员少、证据收集难等困难，给河湖采砂动态监管带来了挑战，急需新的监管手段[4]。

遥感技术具有速度快、成本低、范围广的特点，近年来广泛应用于水体信息提取[6]、水质监测[7]、采砂监管[8-13]等水资源调查、监察、保护领域。目前，遥感技术应用于采砂监管主要包括采砂船识别和高浑浊水体识别2个方面。应用高分辨率卫星[3]或无人机遥感图像[4]，采用人工解译[3-4]或深度学习图像自动识别[8]方法可识别采砂船、采砂场。另外，采砂扰动之后，水体浊度、透明度、总悬浮物浓度、悬浮泥沙浓度等相关指标异常升高，能够通过水色遥感技术判断采砂活动信息[10]。邬国锋等[10-11]建立了鄱阳湖水体透明度(Secchi Disk Depth,SDD)与 LANDSAT TM5 图像红、蓝波段线性经验关系，分析采砂对水体透明度、生物多样性等方面的影响。张伟等[12]利用蓝光、红光、近红外波段计算总悬浮物指数(Total Suspended Matter Index,TSMI)，并建立其与总悬浮物浓度的对数经验关系，分析采砂区时空变化特征。杨静学等[13]应用TSMI分析发现鹤地水库悬浮物浓度高值区存在采砂情况。目前，常规采砂遥感监测需经过大气校正、悬浮物浓度反演两级过程，通过设定总悬浮物浓度经验阈值识别采砂嫌疑区。由于大气及水质组分时空差异大，受限于大气纠正、总悬浮物浓度反演精度，分割阈值常需人工干预，难以实现多时相采砂区自动化识别。此外，总悬浮物浓度反演模型常利用近红外波段泥沙高反射特性，不适用于仅包含红、绿、蓝3个波段的常规无人机自然图像。

水色变化与水中光学组分浓度变化直接相关，通过水色参数可监测水体水质的时空变化，如黑臭水体[14]、水体富营养化[15]、水体透明度[16-17]、水色异常[18]等。主波长水色参数不需近红外波段，仅采用红、绿、蓝3个波段计算，且对大气纠正抗扰性强[17]，在监测水体浊度、河道采砂方面极具应用潜力，而目前这方面的研究鲜见。

本文以枫树坝水库、鹤地水库水域为研究区域，应用Landsat8、Sentinel-2卫星遥感图像开展采砂区主波长水色参数分析，并采用自然断点法自动设定划分阈值，旨在建立一种适用于红、绿、蓝3波段图像，无需大气纠正，自动识别采砂嫌疑区的方法，为常态化采砂监管提供动态监测信息。

2 数据与方法

2.1 研究区介绍

枫树坝水库、鹤地水库均为广东大(I)型水库，因位于省际交界区域，交通不便，管理薄弱，造成非法采砂屡禁不止，难以监管。枫树坝水库集雨面积5 150 km2，总库容19.4亿 m3，位于东江上游的龙川县枫树坝镇，是集防洪、发电、供水功能于一体的综合性水库。近年来，枫树坝水库入库河流寻邬水常出现水体异常浑浊的情况，疑似存在非法采砂现象。鹤地水库集雨面积1 495 km2，总库容11.44亿 m3，库区地跨广东廉江、化州和广西陆川、博白4个县(市)。鹤地水库是一座以灌溉为主，兼顾防洪、发电、旅游的水利枢纽骨干工程，承担着雷州半岛97 333 hm2农田灌溉和近400万人口生活用水供水任务。2012—2018年，鹤地水库交界水域非法采砂行为时有发生，仅依靠人工巡查，难以掌握采砂的时空分布信息。引入遥感监测方法后，通过加大跨区执法力度，并经过长效治理，自2019年之后，鹤地水库非法采砂现象得到了有效控制。

2.2 数据源

采用Sentinel-2卫星T50RLN标准分幅数据，包括2019-11-08、2019-11-23、2019-12-08、2020-02-21、2020-04-16共5个时相图像，于https:∥scihub.copernicus.eu/dhus/#/home网站下载获得，用于枫树坝水库采砂监测。Sentinel-2卫星图像为L1C级产品，是经过几何纠正的表观反射率数据。同时采用Sentinel-2卫星T49QDE标准分幅数据，以及相近时相行列号为124045的Landsat8卫星数据监测鹤地水库，检验方法对不同卫星数据的适用性。Landsat8卫星图像于https:∥earthexplorer.usgs.gov/网站下载，为L1级数据，需进行几何纠正和辐射纠正处理获得表观反射率数据。

2.3 研究方法

2.3.1 CIE-XYZ颜色系统介绍

水体中非藻类悬浮物、浮游植物和有色可溶性有机物，这3种光活性物质及水分子共同决定了水体颜色[17]。由于水中物质组分的物理化学性质不同、浓度不同，水体显示出不同的颜色，因此水体颜色具有表征水体水质的潜力。

为了建立水色与水质的关系，需要对颜色进行数字化表示。国际照明委员会(CIE)于1931年建立了一套标准颜色系统CIE-XYZ颜色系统[16]。该系统将CIE-RGB系统中的3个原色R、G、B转换为3个理想的原色X、Y、Z，关系式[19]为

(1)

CIE-XYZ颜色系统中光谱的三刺激值X、Y、Z全为正值，并且X、Z两原色只代表色度，没有亮度，光亮度只与三刺激值Y成比例[16]。为了便于目视三刺激值定义的颜色，国际照明委员会规定了CIE-xy二维色度图，定义x+y+z=1，对z进行降维处理，使每种颜色都对应一个色度坐标。根据CIE-XYZ系统中的三刺激值X、Y、Z计算色度图上二维坐标(x,y)，计算公式为

(2)

图1 CIE-xy色度图Fig.1 CIE-xy chromaticity diagram

CIE-xy二维色度图如图1所示。图1中S称作等能白光点，其色度坐标为(1/3，1/3)，表示3种原色等量混合。C点是任一颜色色度坐标点，坐标为(x,y)，从S指向C的直线与光谱轨迹的交点D所对应的波长，即该颜色的主波长λd。补色波长λc为从C指向S的直线与光谱轨迹的交点M所对应的波长。从等能白光点到色度图外围线的每一条射线上的点都具有相同的颜色主波长。颜色主波长λd是颜色量化的重要指标，它将可见光颜色从380～700 nm以1 nm间隔来标识，能够以具体的波长形式表示一种颜色的色调。

2.3.2 主波长水色参数采砂识别算法

水库水域在采砂扰动之后，水体浊度升高，水体光谱反射峰向红光通道偏移，水体表现为黄棕色。水体逐渐浑浊，水色定量表征体现为主波长水色参数由蓝绿波段向黄红波段逐渐递增，即水体浊度升高，主波长增加。

原始遥感图像经过预处理后获得表观反射率图像，再依据改进归一化差异水体指数(Modified Normalized Difference Water Index,MNDWI)进行水陆分离[20]，获得水体掩模。将水体掩模范围内红、绿、蓝3波段表观反射率代入式(1)，计算X、Y、Z三刺激值。根据式(2)，将XYZ转为色度图坐标(x,y)。采用张晨朝[21]提出的分区法，由色度坐标(x,y)计算获得水体主波长水色参数图像。

进而采用自然断点分类法Natural Breaks(Jenks)[22]对主波长图像进行分类测试，分析确定高浊度水体的划分阈值，可实现采砂嫌疑区识别。自然断点分类法的思想是在确定分类个数的情况下，计算各种阈值组合下分类的方差和，其值最小的就是最优的分类结果，即可自动确定最优分类阈值。根据已知试验区的情况，主波长图像分为高、中、低3个等级，数值最高的等级可识别为采砂嫌疑区。主波长水色参数计算及自然断点分类程序算法均采用python语言，调用gdal开源库编程实现。

3 结果与分析

3.1 鹤地水库

采用本方法对Landsat8卫星、Sentinel-2卫星多时相遥感图像进行处理，测试不同传感器、不同大气状况下，自动识别鹤地水库采砂嫌疑区的效果，结果如图2、图3所示。

图3 鹤地水库Sentinel-2卫星遥感采砂嫌疑区识别Fig.3 Identification of suspected sand mining area by Sentinel-2 satellite remote sensing in Hedi Reservoir

Landsat8卫星不同时相主波长水色参数图像显示了水体浊度差异，九洲江水域及化州库湾异常浑浊区域主波长值域范围为525～590 nm，浑浊程度较小的库心区域主波长值域范围为480～490 nm。根据自然断点法分别以最高类阈值>525.45、525.50、525.47、525.47 nm将2015年1月1日、2016年11月3日、2017年12月2日、2017年12月24日共4个不同时相主波长高值区识别为采砂嫌疑区。由4个不同时相的结果可知，采砂区划分阈值十分稳定，即使2017年1月22日及2017年12月24日存在薄云，阈值差异仍然非常微小，说明方法抗扰性强，无需对图像进行大气纠正预处理，易于使用。

Sentinel-2卫星与相近时相Landsat8卫星的主波长水色参数图像反映的水体浊度一致，2016年10月30日、2017年1月28日、2017年10月30日、2017年12月19日共4个不同时相采砂嫌疑区识别阈值分别是522.55、522.64、522.47、525.72 nm，差异也非常小。由此说明主波长水色参数对于不同传感器具有很强的适应性和稳定性。

Landsat8卫星2015年和2016年2个时相、Sentinel-2卫星2016年一个时相主波长水色参数遥感结果显示，广西、广东交界九洲江水域异常浑浊，2017年Landsat8、Sentinel-2卫星不同时相的遥感结果均显示，除了九洲江水域，茂名化州与湛江廉江交界库湾水域也异常浑浊。主波长水色参数反映的鹤地水库水域浑浊程度时空分布特征与杨静学等[13]的研究结果相似，均体现了九洲江及化州库湾存在采砂嫌疑。但由于主波长水色参数模拟人眼对水色的感官，将红、绿、蓝3个波段信息压缩，受大气状况影响小，可通过自然断点分类法自动划分高值区识并将其别采砂嫌疑区，具有人工干预少、抗扰性强的优势，更适用于多时相遥感反演结果对比分析。

图4 枫树坝水库Sentinel-2卫星遥感采砂嫌疑区识别Fig.4 Identification of suspected sand mining area by Sentinel-2 satellite remote sensing in Fengshuba Reservoir

3.2 枫树坝水库

对2019年11月8日至2020年4月16日共5个时相Sentinel-2卫星遥感图像计算主波长水色参数，结果如图4所示。采用自然断点法划分主波长水色参数为3个等级，经测试各时相最高值分类阈值非常稳定，为525.4±0.1 nm。以主波长高值区识别为采砂嫌疑区，由图4可见，方法自动性强、稳定可靠，既能识别水体较为浑浊时的情况，也未对2020年2月21日及2019年11月8日水体较为清澈时的情况产生误判。

多时相主波长水色参数遥感结果显示，广东、广西交界寻邬水河段经常出现异常浑浊现象，与周围支流水色形成明显对比，极可能存在采砂现象。根据2019年11月6日以及2020年12月17日高分辨率卫星遥感图像显示，该河段在图4(a)中A、B点位存在采砂船。如图5所示，在B点位清晰可见采砂船周围水体浑浊，正在进行采砂活动。由高分辨率卫星遥感图像采砂船解译结果，证实了主波长水色参数采砂区识别结论。

图5 枫树坝水库采砂船识别验证Fig.5 Verification and identification of sand dredgerin Fengshuba reservoir

4 结 论

经枫树坝水库、鹤地水库不同地区，Landsat8、Sentinel-2卫星不同遥感图像试验发现，主波长能够反映水色，其值随水体浊度升高而增大。采用自然断点法划分主波长为高、中、低3个等级，高值区分类阈值非常稳定。采砂区高浑浊水体可采用自然断点法自动选择阈值判定，或者依据本文试验采用主波长>525 nm作为判断条件。试验证明主波长水色参数在识别采砂方面具有方便快捷、自动识别、适用多种卫星图像及无人机自然图像的优势。方法可为采砂动态监管提供时空分布信息，辅助管理人员有针对性地开展重点巡查，有效提升常态化采砂监管的能力。

由于缺少实测数据，本文尚未建立浊度与主波长水色参数的定量关系，后续需进行光谱采样试验，选取水色参数和光谱指数建立深度学习水体浊度遥感反演模型。此外，可进一步提升数据处理效率，应用Google Earth Engine遥感云平台高频率计算Landsat8、Sentinel-2等卫星图像主波长水色参数，统计多时相采砂嫌疑区识别频度图像，快速评估经常性浑浊水体，以此排除因雨洪影响而浑浊的水域，识别异常情况，提供采砂监管辅助信息。