湖泊是区域生态环境中的重要一环，对流域内居民生命安全和社会经济发展有着巨大的影响. 随着我国经济的发展，湖泊正在受到越来越多的人为因素的干扰，水域萎缩和水质恶化等问题日渐严重，这些都成为了区域可持续发展的瓶颈^[1-4]，采砂便是其中影响较大的人为活动之一. 作为建筑砂石的重要来源，采砂有着较大的经济效益，但也会对水域生态环境产生较大的破坏. 过量的采砂会导致大量的泥沙搅动，改变河道及水底地形，进而影响水下环境及相关动植物生物量^[5-10]. 此外，采砂活动会导致地质疏松，很容易引发泥沙的再悬浮，进而改变水体浊度. 因此，采砂活动对水环境及水生态有着巨大的影响^[11-12]. 我国多个湖泊的生态环境都遭受了采砂活动的影响. 如洪泽湖在2012年后由于过度的采砂活动导致水体悬浮颗粒物浓度增高，生态环境遭到一定程度的破坏^[13]，更有研究表明，洪泽湖采砂区域的营养盐浓度显著高于水生植被生长的区域，透明度和生物多样性却显著较低^[14]. 鄱阳湖作为目前我国第一大淡水湖泊，承纳来自赣江、抚河、信江、饶河和修水汇入的大量泥沙. 统计数据显示，鄱阳湖每年淤积的泥沙量可达1000万t以上，拥有丰富的砂石资源. 受长江中下游河道禁止采砂影响，鄱阳湖丰富的湖砂资源成为了开发的方向^[15-19]，因此，对采砂活动进行监测与管控对鄱阳湖流域的可持续发展有着重要意义.

遥感在采砂活动监测领域开展的应用多与水体悬浮泥沙浓度相结合，典型的研究结果以洪泽湖和鄱阳湖2个湖泊为代表. Cao等利用MODIS影像与实地巡查数据建立了洪泽湖悬浮颗粒物浓度估算算法，并将2012年起洪泽湖悬浮颗粒物浓度的显著变化归因于密集的采砂活动^[13]. Duan等则在日间影像的基础上使用了VIIRS昼夜波段(DNB)夜间灯光(NTL)数据，总结得到了洪泽湖2012-2017年采砂船的时空分布特征，并发现洪泽湖悬浮物浓度呈现的变化规律与采砂船的时间分布相一致^[14]. 针对鄱阳湖采砂的遥感研究相对更多一些，如de Leeuw等利用卫星图像估计了禁砂期离开鄱阳湖的船只数量，并评估了鄱阳湖采砂的影响^[20]，以及Lai等发现鄱阳湖的外流河道会随着采砂强度的变大而加宽加深^[21]. 较为全面系统的研究结果以Li等和Feng等^{[15, 21]}为代表. Feng认为，大量涌入的采砂船是影响鄱阳湖水体悬浮泥沙浓度的一大驱动因素，并建立了鄱阳湖北湖年平均悬浮泥沙浓度与北湖最南端两段船舶数的回归关系^[16]. Li通过湖面分区、船舶目标增强和悬浮泥沙反演结果叠加分析来提取鄱阳湖采砂船作业点，并用底质类型分类结果进行结果的验证^[15]. 不过，该研究结果覆盖的时间较早，对2013年以后鄱阳湖采砂活动的变化规律未见较为系统的研究结果. 本文使用分辨率更高的GF-1 WFV数据，结合采砂船常与运砂船并列排布以及对附近水体透明度有巨大扰动这两大特点^[22-23]，对2013-2020年鄱阳湖卫星影像进行船舶增强处理和水体透明度反演，识别鄱阳湖上的采砂船作业点，进而分析鄱阳湖采砂活动的时空间分布规律及其驱动因素.

1 研究区域与数据源 1.1 研究区概况

鄱阳湖(28°22′~29°45′N，115°47′~116°45′E)(图 1)地处我国江西省北部，是目前我国水域面积最大的淡水湖泊. 鄱阳湖区域属亚热带湿润季风型气候，冬春季气温较低，夏季潮湿多雨，秋季晴热干旱^[24]. 年平均气温在17℃左右，年内变化较大. 年日照量约为1800~2100 h，雨量丰沛，年平均降水量可达1400~1900 mm，但降水年际变化较大^[25-26]. 鄱阳湖的湖区一般以都昌县南的松门山岛为界分为南北两个部分^{[15, 27]}. 北部湖区狭长且直通长江，为入江水道；南部湖区在丰、枯水期差异较大，丰水期湖面辽阔，枯水期水道狭长，具有独特的“洪水一片，枯水一线”特点. 整个湖区的水量呈季节性变化，一般来说，丰水期为4-9月，湖水淹没草洲、滩涂，主湖区面积辽阔；枯水期为10月至次年3月，水位较低，草洲、滩涂重新显现. 鄱阳湖是采砂船良好的作业区域^[28-29]，其采砂活动最早见于1990s，2001年开始兴起，乱采滥挖现象日渐严重^[15]. 2006年起，相关部门陆续出台政策对采砂活动加以限制. 江西省水利厅(http://slt.jiangxi.gov.cn/)和九江市采砂管理局(http://csglj.jiujiang.gov.cn/)发布的政策文件显示，在2008年全面禁止采砂一年后，鄱阳湖自2009年起采取批复式采砂管理办法，根据每年情况确定可开采量.

1.2 研究数据

本文选用高分辨率的GF-1多光谱相机(WFV)数据对鄱阳湖采砂船作业点进行识别. GF-1卫星于2013年4月开始提供数据服务，WFV数据的空间分辨率为16 m，重访周期2 d，幅宽800 km，搭载红、绿、蓝光及近红外4个波段，波长范围分别为0.45~0.52、0.52~0.59、0.63~0.69和0.77~0.89 μm. 由于鄱阳湖承担了流域内较多的航运功能，河道上常年布满各色的渔船和运输船，使用空间分辨率更高的卫星影像可以更加清楚地对船只进行识别，减轻验证环节的压力. 此外，GF-1 WFV在可见光波段的设置与HJ-1A/B卫星传感器的波段设置基本相同，便于本研究直接选用发展成熟的模型. GF-1 WFV数据的不足在于幅宽较窄，一景影像有时并不能覆盖鄱阳湖湖区的全境. 但现有研究显示，鄱阳湖采砂活动绝大部分集中在28°55′~29°35′N范围内，一景GF-1 WFV影像的幅宽能够满足这一要求. 因此，本文使用GF-1 WFV影像对鄱阳湖采砂船作业点进行监测.

本研究共下载了149景2013年5月-2020年12月鄱阳湖区影像，根据云层覆盖位置和云量进行了筛选，共余下133景有效数据. 使用ENVI5.3中的RPC Orthorectification Using Reference Image工具对其进行几何校正，参考相近成像时间的Landsat 8 OLI影像进行参考点的自动采集与校正. 大气校正方面，使用FLAASH Atmosphere Correction工具消除Ⅱ类水体中水汽和气溶胶散布的影响. 鄱阳湖所处的纬度范围在30°N附近，因此大气模型在7-9月选择Tropical，其余时间选择Middle-Latitude Summer. 气溶胶模型选择Urban，Initial Visibility采用默认的40 km.

2 研究方法

经团队调研发现，鄱阳湖湖区的采砂船主要为大型抽砂泵式采砂船，长70~80 m，宽15~20 m. 由于采砂船本身不具备运载砂石的能力，因此常与多艘运砂船一同作业，呈并列排布. 尤其在丰水期时，辽阔的湖面甚至可能有数艘采砂船与十数艘运砂船组成的巨大采砂平台，在GF-1 WFV影像上约为10×6个像素左右. 同时，根据实地考察和Google Earth高清影像的结果可以看出，这些船舶多为许多颜色组成，与运输船仓底露出的单一色调有明显的差异，在影像上具有斑驳的纹理特征. 这从形状和纹理特征的角度提供了区分采砂船作业点与湖区内其他渔船、运输船的方法. 典型采砂船在16 m分辨率的GF-1 WFV影像上如图 2所示.

同时，正在作业的抽砂泵会造成附近水体的扰动，使沉积的泥砂再次悬浮，进而导致附近水体的悬浮泥沙浓度异常增高. 鄱阳湖水生植被较少、蓝藻水华暴发的频率也很低，悬浮泥沙浓度的升高是导致水体透明度下降的主导因素. 因此，本文判别采砂船的主要依据有二：一是通过船体的大小及排布形态筛选初步的采砂船疑似结果，二是在此基础上通过提取水体透明度的突变区域判断其是否在工作. 此外，采砂船对作业区域的湖底底质类型有着较严格的要求，即必须以沙滩或泥滩为主^[30]，且较少出现在鄱阳湖航船的航线上. 这些也是本文判别采砂船的辅助手段.

2.1 船舶目标影像增强处理

一般来说，船舶目标与水体的反射率在近红外波段(GF-1 WFV第4波段)的差异最大，但鄱阳湖水体含泥沙量较高，仅用近红外波段区分高浑水体与船舶目标的效果易受到影响. 考虑到船舶目标在绿光波段到近红外波段的反射率变化趋势相较于水体要平缓很多，本文在近红外波段的基础上引入绿光波段，使用归一化差值水体指数(normalized difference water index，NDWI)来突出船舶目标^[31]，公式如下：

$ N D W I=\frac{R_{\text {green }}-R_{\mathrm{NIR}}}{R_{\text {green }}+R_{\mathrm{NIR}}} $

(1)

式中，R_green为绿光波段反射率，R_NIR为近红外波段反射率. 这一算法可以增强船舶目标与水体的差异，还可以用于后续的水体提取. GF-1 WFV原始影像及波段比值增强后影像效果对比如图 3所示.

2.2 水体透明度反演

水体透明度能够反映水体能见度和光学性质，与水体组成成分有较高的相关性，利用遥感影像反演水体透明度的研究已经屡见不鲜. 本研究使用Li等结合HJ-1A/B CCD遥感影像与实测透明度数据建立的适用于鄱阳湖的水体透明度(secchi disk depth，SDD)反演模型^[15]，公式如下：

$ \ln (S D D)=1.454 \frac{R_{560}}{R_{660}}-25.014 R_{830}-2.20 $

(2)

式中，R₅₆₀、R₆₆₀和R₈₃₀分别表示HJ-1A/B CCD在第2、3、4波段的反射率. 该模型在时、空间上都较为稳定. 由于HJ-1A/B CCD与GF-1 WFV的波段设置基本相同，且这一步的目的只需要反演结果能反映水体透明度的突变现象，因此将其应用到本研究中. 该水体透明度反演模型已经应用在Li等针对鄱阳湖采砂活动的研究，效果较好，平均绝对误差(mean absolute error，MAE)为21.6%，均方根误差(root mean square error，RMSE)为30%，R²=0.81.

由于鄱阳湖水域广阔且水质分布不均匀，不适合使用单一阈值对水体透明度的突变区域进行提取. 而Sobel算子可以更好地识别出影像中与周围像元相对差异较大的部分. 因此，本研究将检测出船舶目标的小块均匀水体作为输入，利用Sobel算子提取水体透明度的突变区域，结果如图 4所示.

2.3 其他判别手段

由于鄱阳湖具有“洪水一片，枯水一线”的特点，大面积的草洲、洲滩等在枯水期露出水面，可以利用这一特点对鄱阳湖的底质进行监测与分类. 由于湖底底质类型的变化较慢，因此这里参考了刘子潇的底质类型分类结果图^[30]. 采砂作业常在沙滩区域进行，泥滩区域偶有出现，草洲/芦苇区域不适合进行采砂作业，这些区域上应避免出现采砂船的提取结果.

此外，鄱阳湖承担了流域内较多的航运任务，湖面上常有各类航船游弋. 而采砂船在短时间内的工作范围变化相对较小，因此一般会避开航线以免相互妨碍. 鄱阳湖上航船的航线相对固定，在多景影像上有较好的体现，依据2013年10月5日的影像提取出的大致航线如图 5所示. 在提取采砂船时，应尽量避免处在航线上的结果.

3 结果与分析 3.1 总体结果

在对现有的133景鄱阳湖GF-1 WFV影像进行处理后，在其中的117景上检测出了采砂船作业点，且绝大部分集中在主航道两侧. 对比鄱阳湖底质类型分类结果图可以看出，识别出的采砂船与底质分析结果中最适合采砂船作业的区域——沙滩相对应，还有部分采砂船作业点在底质为泥滩的区域被识别出来. 而被认为不适宜采砂船作业的草洲和芦苇区域(大部分分布在湖区南部区域)则基本没有识别出采砂船. 识别出采砂船的区域与团队实地考察的结果也相符，这说明本研究的结果整体上是可信的. 限于篇幅，这里仅选取了3景影像及其识别结果(分别代表枯水期、丰水期和平水期)进行了展示，如图 6所示.

2013-2020年鄱阳湖采砂船统计结果见表 1，需要说明的是，鄱阳湖的采砂船数量在2013-2018年仍然呈现出明显的分期特点：即丰水期和枯水期提取的采砂船数呈现出较大的差异. 因此在之后列出的采砂船统计结果中，按丰水期和枯水期加以区分. 2019年后，采砂船检测数在丰/枯水期的差异不再明显，丰/枯水期采砂船的数量相较之前显著减少.

3.2 采砂船提取结果验证

为确保本研究采砂船提取方法的可行性，对上述提取的采砂船位置进行了验证. 2019年10月2日星子县附近以南水域上GF-1 WFV影像及采砂船识别结果如图 7a所示，验证数据为2019年10月4日Google Earth高清影像，船只的轮廓与纹理清晰可见，如图 7b所示.

图 7中，(a)图的红色实心三角形标记为2019年10月2日GF-1 WFV影像提取的采砂船作业点，(b)图的空心三角形标记为2019年10月4日Google Earth高清验证影像识别出的采砂船，忽略2天时间间隔带来的采砂船位置偏移. 对比两次结果，GF-1 WFV影像提取出4个采砂船作业点，Google Earth影像则识别出6个. 2个漏提的采砂船作业点(图中黄色空心三角形标记)分别出现在左下角和中部. 左下角漏提采砂船是在后期总结时被误归至邻近采砂船作业点中，将2个作业点识别为1个. 中部采砂船的漏提是由于水体透明度的突变过程没有得到很好的反映，可能是受到附近船只较多的影响，同时也不能排除采砂船未在作业的情况. 可以看出，基于经验的人为提取过程和附近过于密集的船舶排布会对采砂船作业点提取的准确度带来一定的影响. 但其余4个提取出的采砂船作业点(图中红色空心三角形标记)与验证影像的结果吻合度很高. 总的来说，除了船舶数量过于密集的区域可能出现部分漏提的情况外，大部分提取结果是准确的，可以认为本研究的采砂船提取结果能正确地反映采砂船的数量及分布情况.

3.3 鄱阳湖采砂作业点时间变化规律

2013-2020年间鄱阳湖的采砂活动强度呈先增后减的趋势. 2013-2016年平均每年检测出的采砂船数逐年增多，2016年后开始减少，2019-2020年再次锐减. 近年来，鄱阳湖水体面积萎缩、湖区干涸等水环境问题时有发生，大规模的采砂作业及其带来的环境变化已经引起了当地政府部门的注意，政策调控随之成为鄱阳湖采砂活动最大的影响因素^[13-15]. 查阅相关资料发现，2008年起江西省水利厅对鄱阳湖地区的采砂活动采取了科学考量、批复开采的政策，每年的采砂活动都要受到当年可开采量的限制. 查阅江西省水利厅公布的鄱阳湖2013-2018年批复可采砂量时发现，这几年采砂船检测数的年际变化与当年批复的可采砂量有着很高的相关关系，两者的变化趋势吻合度极高，如图 8所示. 在采砂船检测数明显较多的2016年，政府批复的可采砂量也明显多于其他年份. 而随后的2017年和2018年，检测到的采砂船数也随着批复采砂量的大幅降低明显减少. 2019年起，鄱阳湖采砂政策转变为更加严格的全面禁止采砂，影像上检测到的采砂船数也随之锐减. 因此可以认为，政策调控是2013-2020年间鄱阳湖采砂船数量变化最主要的影响因素.

在采砂船检测数的年内变化方面，由于月份易受到可用影像数、当年汛期时段等因素的影响，本研究使用各年采砂船提取数最值出现的季节作为指标说明规律. 结果显示，2013-2020年鄱阳湖每年采砂船数量的高峰期均出现在夏季，秋季也偶有出现，低谷期均出现在冬季，春季也偶有出现. 总的来说，采砂船数的年内变化呈现出先增后减的趋势，如图 9所示，即春季起采砂活动逐渐活跃，在夏季达到顶峰，之后到冬季一直下降. 值得一提的是，2019-2020年的采砂船检测数虽然仍呈现夏季高冬季低的规律，但年内变化明显地趋于平缓，极值也更多地出现在春、秋两季.

3.4 鄱阳湖采砂作业点空间变化规律

采砂作业会搅动水中的悬浮泥沙，对一定范围内水体的水质和透明度都有较大影响. 由于这一影响效果难以量化，本研究使用距离来说明采砂作业点对附近水体的影响，将采砂船作业点识别结果的点状数据转化为面状数据. 经过之前的实地调研发现，鄱阳湖上多为抽泵式采砂船，其工作1 d移动的距离约为20 m. 而本研究使用的丰水期影像时间跨度约为25 d，考虑到鄱阳湖丰水期辽阔的水域面积，在理想状况下采砂船连续工作25天的移动距离约为500 m，故而本研究将丰水期采砂船的第一级影响范围设为500 m，第二、三级分别扩大为1000和2000 m. 枯水期水域面积大幅减小，水道狭窄，因此将各级影响范围的分别下调为300、600和1000 m. 结果如图 10所示.

丰水期识别到的采砂船主要分布在星子县、笔架山附近和松门山岛附近水域，也有部分采砂船在笔架山至松门山岛中间的水域和周溪镇附近被识别出来，在北部的屏峰村水域也有零星的分布. 而在枯水期，采砂船主要分布在松门山岛附近和星子县、笔架山附近水域. 总的来说，鄱阳湖的采砂活动主要集中在松门山岛以北、周溪镇以及星子县、笔架山附近水域. 而从前面的研究结果可以看出，2019-2020年，鄱阳湖采砂船作业点的数量在丰/枯水期的差别几乎消失，在位置分布上呈现出相同的特点. 因此，这两年的结果分析中不再进行丰/枯水期的区分(图 11). 识别出的采砂作业点分布在松门山岛附近、星子县以南及笔架山以西水域，但周溪镇附近未见采砂船作业点，南湖区域的采砂活动基本销声匿迹.

此外，Li等的研究指出，在2008年开始采取采砂管控措施后，鄱阳湖采砂活动有向南湖扩展的趋势^[15]. 为探究2013-2020年这一规律是否仍在延续，进一步获取鄱阳湖采砂活动的空间动态，本研究将检测到的所有采砂船纬度进行了统计，如图 12所示. 可以看出，鄱阳湖采砂活动在纬度分布上明显向北湖集中. 2013-2017年采砂船出现的纬度范围较大，南部湖区检测出的采砂船数量较多，2018年后采砂船出现的纬度集中在北部湖区. 对各年在南湖检测出的采砂船平均数进行统计后发现，其呈现出先增后减的趋势，在2018年后控制在很低的水平. 这些都表明，鄱阳湖采砂活动向南部湖区扩张的趋势已得到遏制.

4 结论与展望

本文利用采砂船工作时会明显改变附近水体的悬浮泥沙浓度这一特点，结合采砂船的形状和纹理特征、鄱阳湖底质类型及航线位置对2013-2020年鄱阳湖采砂活动进行了监测，并总结分析了鄱阳湖采砂活动的时空间变化规律. 结果显示，2013年以来，鄱阳湖的年平均采砂船检测数呈先增多后减少的趋势，2016年是采砂活动的顶峰，这与鄱阳湖每年批复的可开采量呈现基本相同的变化趋势，2019年后采砂船检测数的锐减也与当年全面禁止采砂的政策相对应. 这说明监管政策是影响鄱阳湖采砂活动最主要的因素. 同时，鄱阳湖采砂活动向南湖扩展的趋势也得到了有效的遏制. 不过，在2019年全面禁止采砂后，仍然在松门山岛、星子县及笔架山附近监测到了零星的采砂船作业点. 本研究可以评估近年来鄱阳湖采砂活动监管的成效，并给相似水域采砂活动的监管及其效果评估提供参考.

本研究也能为计算机智能监测采砂活动提供实现流程. 虽然本研究在人工筛选采砂船结果的过程中会产生一些误差，在船只密集的区域有时会出现漏判、错判的情况，且最终结果易受到影像分辨率、船舶阴影及云雾遮挡的干扰，但是随着更多性能优异的卫星传感器的发射、影像处理算法的不断进步和机器学习应用的日渐成熟，更为精确、便捷的采砂船计算机处理方法有望在将来得到开发和应用.