基于遥感技术对松花江悬浮泥沙浓度的分析

韩 旭，杜 崇，刘福全，李 瑞

(黑龙江大学 水利电力学院，哈尔滨 150080)

0 引 言

河流的悬浮泥沙是河流重要的水质参数之一，泥沙输移是河流中重要的水文现象。水体悬浮泥沙浓度影响一系列的水体光学特性，如河流水色、水体浑浊度、水体透明度等，进而影响水生生物的生态环境，例如光合作用、初级生产力、营养流动、河流生物多样性等。在天然江河修建水坝，水坝上游水位抬高，流速变缓，河流中的悬浮泥沙在坝底淤积，泥沙压力可能造成水坝应力分布不均，破坏坝体结构和稳定性。因此在泥沙含量偏高的河流中修建水坝，应考虑到应力的影响。此外，河流泥沙含量对河道平衡也有很大影响，局部泥沙淤积造成河床变形等问题。因此对河流泥沙的动态监测有重要意义[1]。传统测量悬浮泥沙的方法是对采样区内的采样点进行逐个点的长时间监测，这种方法不仅耗时过长，而且存在着成本高和无法对大面积水域进行长时间监测等多方面缺点。遥感技术的监测范围广、耗时短，且遥感图像获取的成本较低，运用遥感技术对河流悬浮泥沙进行监测成为新的发展趋势。受制于遥感图像分辨率等多方面因素，遥感技术对悬浮泥沙的反演多是对海洋和湖泊等面状水体，对河流线状水体的研究较少。韩震等[2]对遥感卫星单波段的反射率和实测悬浮泥沙浓度构建了数学模型对其进行研究，但单波段的反演模型无法排除水中叶绿素等物质的干扰。方馨蕊等[3]使用随机森林回归模型对悬浮泥沙浓度进行遥感估算，但机器模型需要大量的训练样本，很难在河流遥感的定量反演中应用。林承达等[4]利用LandsatETM+遥感影像对长江中游悬浮泥沙浓度进行遥感反演，将遥感技术运用到线状水体悬浮泥沙的监测上，但LandsatETM+影像分辨率为30 m，对于线状水体，其分辨率较低。本文使用Sentinel-2卫星数据，其遥感影像分辨率经过SNAP软件重采样后分辨率为10 m，且相较于Landsat系列陆地卫星，Sentinel-2卫星重访周期更短[5]。

1 数据与方法

1.1 研究区域概况

由于大顶子山航电枢纽工程，大坝上游松花江哈尔滨段水位已蓄至115 m，两岸大面积的平地和滩地被回水覆盖，形成从枢纽大坝至肇源县长达128 km的人工湖。大顶子山航电枢纽工程蓄水后，松花江左岸原有的湿地资源被成倍放大。位于呼兰区的“黑龙江呼兰河湿地自然保护区”申请立项，经批准的保护区总面积为19 262 hm2，其中林地1 680 hm2，芦苇湿地、湿生草地、草甸等近万hm2，是国内目前最大的城区湿地。大顶子山航电枢纽建成后，松花江哈尔滨段的通航期由原来的100 d左右提高至210 d左右，在枯水期时，1 000 t船队也可以在松花江哈尔滨段畅通无阻。对于河流遥感的定量反演，河面宽度增大，卫星传感器足够清晰识别河面目标，适合作为本次实地采样区域。

近年来，松花江大力发展航运，大顶子山枢纽附近的水运量在1×107t以上，研究该区域的泥沙分布，不仅为航道建设、航路规划提供数据支持，也为防止大顶子山航电枢纽泥沙淤积提供参考。

1.2 实测数据

枯水期和丰水期野外实地采样时间分别为5月12日，7月13日。两次采集的采样点均是14个，采样点选取时，选择江面平缓的位置，并且当日无降水，反演模型建立点均匀分布在河流横断面上，确保反演模型的准确性。14个采样点中，1～8号点用于反演模型的建立，9～14号点用于模型精度的验证。使用手持GPS进行定位，确保两次采集点位置相同，悬浮泥沙采样方法与实验室测量方法参照国家河流悬移质泥沙监测规范(GB/T 50159-2015)。采样点位置见图1。

图1 采样点位置

实地测量前先测量每个采样点对应容器和滤膜的质量，采用容器为500 mL宽口塑料瓶，因为各个波段对水体的穿透性有限，现场采样水深为0.5 m，每个采样点采集3瓶水样取平均值，目的是为了减少误差和后期数据筛选。采集水样后，即到实验室进行处理。枯水期和丰水期各采样点悬浮泥沙浓度见表1和表2。

表1 枯水期实测采样点悬浮泥沙浓度

表2 丰水期实测采样点悬浮泥沙浓度

1.3 遥感数据

本文使用的遥感数据均来自Sentinel-2卫星，该卫星遥感数据在欧空局网站(https:∥scihub.copernicus.eu)下载，Sentinel-2卫星于2015年6月23日升空，遥感图像宽幅为290 km，重访周期为5 d，拥有13个波段，各波段参数见表3。

表3 Sentinel-2卫星各波段参数

实测时间为5月12日下午14：00左右，7月13日上午10：00左右。由于Sentinel-2卫星图像标注时间为格林尼治时间，格林尼治时间与北京时间相差8 h，所以5月12日的Sentinel-2卫星图像时间选择尽可能接近上午6：00，7月13日的Sentinel-2卫星图像时间选择尽可能接近凌晨2：00，减小由时间差异造成的误差。最终选择5月12日上午6：32的图像和7月13日上午4：17的图像，该图像云层覆盖率分别为3.12%和9.22%，均小于10%，能见度高，适合运用于反演模型的建立。

因Sentinel-2卫星经过处理的L2A数据没有与本次实际采样时间对应的遥感影像，本次Sentinel-2卫星数据下载的是L1C级产品，经过几何校正，没有进行辐射定标和大气校正，需要进行遥感数据的预处理。大气校正方法包括ENVI中FLASSH模块大气校正、6S模型大气校正和Sen2cor等多种方法。由于本文选用Sentinel-2卫星遥感图像，选用由欧空局提供的Sen2cor将L1C级数据进行处理，得到L2A级大气底层反射率数据[6]。

Sentinel-2卫星各个波段的空间分辨率不同，需要进行重采样将空间分辨率统一，运用SNAP软件将各波段的空间分辨率重采样为10 m，重采样方法为双线性插值法。输出数据为ENVI可以打开的格式，运用ENVI中的LayerStacking工具进行波段融合后，进行图像裁剪和拼接处理，获得各个波段实测点的遥感反射率数据。

2 反演模型构建

2.1 敏感波段选取

Senlinel-2A拥有13个波段，第10波段用于大气校正，不参与相关性研究，从第1波断开始，中心波长不断增加，第1波段到第12波段其波长为0.433～2.190 μm，由于松花江相较于国内其他水系是少沙河流，在较清澈水体中，河流的反射率和波长大致呈反比关系，即波长越大，反射率越低。通过影像分析和波谱分析，在波长1.0～2.190 μm，泥沙较少的自然清澈水体反射率极低，近乎为0，可以忽略不计。此外，Band5、Band6、Band7是植被红边波段，对水体的敏感性较低，最终选择Band2蓝光波段、Band3绿光波段、Band4红光波段、Band8近红外波段作为研究波段[7]。

2.2 波段相关性分析

目前用于悬浮泥沙反演的波段组合有单波段、波段比值、多波段组合等多种方法。波段比值和多波段组合，可以减轻水体叶绿素和其他悬浮物质对反演结果的影响。本文利用ENVI中的Band math工具计算单波段、波段比值、3波段组合、4波段组合数据，并与悬浮泥沙浓度做Pearson相关性分析[8]。各波段与悬浮泥沙的相关性系数见表4。

由表4可见，单波段的相关性都比较低，尤其Band8近红外波段的相关性在枯水期和丰水期的相关性系数只有0.096和0.041，相关性极低。只有Band3绿光波段和Band4红光波段的P值小于0.5呈正相关性，并且其相关性系数不大。在波段比值分析中，总体相关性高于单波段相关性，枯水期的B2/B3、B2/B4、B3/B4的相关性大于0.5，丰水期的B2/B3、B2/B4的相关性大于0.5，表现良好。在多波段的相关性分析中，相较于波段比值和单波段其相关性稍好，枯水期的B2/(B3+B4)、B3/(B2+B4)相关性大于0.6，丰水期的B3/(B2+B4)、(B2+B3)/B4相关性大于0.6，可以满足反演模型建立的相关性需求。总体来说，无论单一Band8近红外波段，还是波段比值和多波段组合中涉及Band8近红外波段，其与悬浮泥沙的相关性都比较低。在对黄河悬浮泥沙含量进行监测时，近红外波段是一个不可或缺的波段，但对松花江悬浮泥沙监测时相关性极低。分析原因，可能是由于松花江是少沙河流，Band8近红外波段对其反应不敏感，或是采样区域叶绿素等有机质含量过多，干扰了Band8近红外波段对悬浮泥沙的监测效果[9]。

表4 各波段与悬浮泥沙相关性系数

2.3 反演模型建立

根据波段相关性分析，本文选取相关性系数大于0.6的波段组合建立反演模型，反演模型选用线性模型、指数模型、二次多项式模型、对数模型。运用SPSS曲线估算模块进行反演模型的建立，遥感反射率作为自变量，悬浮泥沙浓度作为因变量[10]。选出拟合度最高的模型作为最终的反演模型。枯水期和丰水期具体模型关系式和拟合度分别见表5和表6。

表5的拟合度分析结果表明B3/(B2+B4)的二次模型拟合度最高，拟合度为0.904，表6的拟合度分析结果表明B3/(B2+B4)的对数模型拟合度最高，拟合度为0.832，基于此结果，提出松花江枯水期和丰水期悬浮泥沙浓度的反演模型分别为

表5 悬浮泥沙浓度经验模型(枯水期)

表6 悬浮泥沙浓度经验模型(丰水期)

f(x)=0.589x2-0.291x+0.124

(1)

f(x)=1.118lnx+0.112

(2)

2.4 模型精度验证

对松花江悬浮泥沙浓度反演模型的精度评价，采用平均相对误差(MAPE)和均方根误差(RMSE)两种方法。

1)平均相对误差(MAPE)：求出每个点相对误差的绝对值后取平均值，平均相对误差越小证明模型精度越高，反之平均相对误差越大则模型精度越低。

(3)

2)均方根误差(RMSE)：又称标准误差，是模拟值与实测值之差的平方与测量次数n的比值的平方根，均方根误差越小表明模型精度越高，反之均方根误差越大表明模型精度越低[11]。

(4)

式中：n为样本数量；yi为实测悬浮泥沙含量；fi为模拟悬浮泥沙含量。

本次模型验证运用实测采集14个点中的6个验证点，在枯水期，6个验证点中最大相对误差为27.17%，最小相对误差为12.3%，平均相对误差为17.42%，均方根误差为0.034 g·L-1。在丰水期，6个验证点中最大相对误差为26.08%，最小相对误差为9.8%，平均相对误差为16.02%，均方根误差为0.016 g·L-1。模型精度基本满足要求。实测值与模拟值的对比见表7。

由表7可见，模拟值总体高于实测值，除个别点外，模拟值的大小与实测值的大小基本呈拟合状态。造成误差的原因：①遥感图像时间与实测时间有近8 h的时差，河流流速等水文条件有所变化[12]；②遥感图像有少量云层覆盖，进而影响后续大气校正结果；③河流中含有叶绿素等物质，影响敏感波段遥感反射率；④实验过程中的过滤，烘干等操作也会造成误差[13]。

表7 验证点模拟值与实测值

2.5 反演效果图

得出枯水期和丰水期的反演模型并进行模型精度验证后，运用ENVI软件制作出枯水期和丰水期悬浮泥沙浓度的反演结果见图2。

图2 反演效果

3 结 论

丰水期的悬浮泥沙浓度总体大于枯水期的悬浮泥沙浓度，其主要原因是丰水期水位上升，将岸边大量的泥沙带入河中，致使河流中的泥沙浓度上升。

悬浮泥沙浓度在河流中的分布总体呈河中心低，岸边高的趋势，其主要原因是岸边流速较低，泥沙容易淤积，且影响遥感反射率的其他杂质较多，反演效果图呈现岸边的悬浮泥沙浓度偏高。