楚森森 ，程亮 , ，程俭，张雪东，刘晋铭
( 1. 南京大学 地理与海洋科学学院，江苏 南京 210023；
2. 江苏省地理信息技术重点实验室，江苏 南京 210023；
3. 中国南海研究协同创新中心，江苏 南京 210023；
4. 自然资源部国土卫星遥感应用重点实验室，江苏 南京 210023；
5. 军事科学院国防工程研究院，北京 100036)
浅海水深数据是航行安全、工程建设、资源开发、海洋救护、生态保护等应用的基础保障[1-3]。传统的船载声呐和机载激光测深方法可获取高精度浅海水深数据，但是成本昂贵，不适用于大规模、周期性测量，并且无法应用于存在争议或被他国非法侵占的岛礁海域[4-6]。
随着卫星遥感技术的发展，利用卫星影像反演水深已成为水深测量的重要手段。卫星遥感有不受地理空间约束、费用低等优势，很好地克服了传统方法的局限性[7]。卫星遥感水深反演的物理基础是可见光具有一定穿透水体的能力，根据比尔-朗伯特定律可知，可见光通过水体时因水柱吸收和散射产生衰减，利用光衰减特性建立光在水中的传输路程（水深）与卫星影像辐亮度关系[8]。卫星遥感水深反演误差普遍为1～3 m[9]，虽然精度不能满足我国海图测量标准，但遥感水深反演方法可解决特殊海域水深信息有无的难题[10]，因此备受关注。目前已有机构实验性地出版基于卫星遥感反演的海图，如英国海道测量局已发布一版基于卫星遥感水深反演的安提瓜岛海图[11]，未来，如何进一步提高卫星遥感水深反演精度仍然是研究重点。
卫星遥感水深反演模型包括解析模型和统计模型。解析模型不需要实测水深数据参与，主要利用水体光传输机理建立反演模型[12]，但模型复杂且未知参数较多，导致反演结果不确定性增加，这些问题阻碍了解析模型的推广应用[13]，简化模型并提高普适性是当前解析模型的研究热点[14-16]。不同于解析模型，统计模型基于少量实测水深点标定模型，建立水深与遥感反射率之间统计关系，因具有结构简洁易操作的优点而被广泛应用[17-18]。统计模型中比较经典的算法有Stumpf比值算法[19]和Lyzenga多项式算法[20]，虽然两种算法参数由统计模型的实测水深样本训练确立，但其数学模型可通过解析模型结合经验公式推导出来[19-21]，因此相对于神经网络、随机森林等算法，Stumpf比值算法和Lyzenga多项式算法，对实测水深样本数量和分布要求较低，且可以一定程度上抑制底质异构、水体异质带来的干扰，算法可靠性与迁移性较强[22]。特别是最近两年，最新研究表明ICESat-2激光卫星数据可以提供少量的实测水深数据用于标定反演模型[23-24]，解决了统计模型算法需要实测水深数据的局限，因此Stumpf比值算法和Lyzenga多项式算法备受关注，并诞生了大量改进算法[25-30]。这些改进算法可归纳为两类，一是非线性拟合改进，即从机器学习角度考虑，如将Stumpf波段对数调节因子由1个增加到2个[31]，将一次多项式模型提升为多次多项式模型[32-36]，将线性拟合改进为支持向量机拟合等[37-38]；
二是地理自适应改进，即从地理学角度考虑，如顾及底质类型水深反演[39]、水深分段反演[40]、地理分区反演[41]等。
统计模型的特点是无需理解水体光传输理论，因此，当前统计模型中经典算法的改进主要从机器学习和地理学角度开展，没有重视光谱测深性能，光谱的测深极限与适用区间会影响水深反演精度，如果进一步从光谱测深性能角度考虑，统计模型反演算法仍有较大精度提升空间。不同光谱测深极限不同，对于清澈水体，蓝光（波长：440～540 nm）最大穿透深度接近30 m、绿光（波长：500～600 nm）最大穿透深度接近15 m、红光（波长：600～700 nm）最大穿透深度接近5 m[42]，但现有水深反演算法没有顾及不同光谱测深特性差异，在不同水深处采用一套反演系数，如在红光可达的浅水区与红光不可达深水区，均赋予红光波段相同的反演系数或权重，影响水深反演精度。对此，本文从光谱测深性能角度改进Stumpf比值算法和Lyzenga多项式算法，提出一种基于多光谱遥感影像的光谱分层策略和一种基于光谱分层的浅海水深遥感反演方法，以我国南沙群岛长线礁和美属维京群岛巴克岛礁为实验区，与经典Stumpf比值算法和Lyzenga多项式算法进行对比分析，证明了本文方法的有效性。
实验区为我国南沙海域长线礁和美属维尔京群岛巴克岛（图1），遥感影像数据为Sentinel-2 L2A级多光谱卫星影像（图1a，图1b），长线礁、巴克岛的Sentinel-2影像获取时间分别为2018年5月9日2时35分（UTC）和2019年10月17日14时57分（UTC），影像坐标系统分别为WGS84 UTM 49N和WGS84 UTM 20N，采用的Sentinel-2波段为近红外、红、绿、蓝波段，空间分辨率为10 m。
实测水深数据包括船载声呐测深数据和机载激光测深数据，其中，长线礁实测水深数据采集仪器为Odom Hydrotrac II单波束声呐设备和SONIC 2024多波束声呐设备，采集时间为2018年1月18-20日，水深基准面为理论深度基准面，水深范围70 m以浅。巴克岛实测水深数据来自美国国家海洋和大气管理局海岸海洋科学中心公开数据集，采集仪器为Laser Airborne Depth Sounder Mark II 激光雷达设备，采集时间为2011年2月21-22日，激光点间距为3 m×3 m，水深基准面为平均低低潮面，水深范围50 m以浅。测深数据垂直精度和水平精度均符合第五版国际海道测量标准1等测量精度，即垂直定位不确定度小于[0.52+（0.013×水深）2]1/2m，水平定位不确定度小于（5+0.05×水深）m。经查潮汐表，长线礁和巴克岛影像获取时刻潮高分别为1.48 m和0.38 m，实验区实测水深已校正为影像获取时刻瞬时海面水深，并提取0～30 m水深作为本文实验水深数据（图1）。水深训练样本和测试样本选用Sentinel-2影像像元均值水深，即对单个像元内部的实测水深数据求均值作为该像元的水深值，长线礁水深像元共计1 127个，巴克岛水深像元共计357 931个。
图1 实验区与水深数据Fig. 1 Study areas and depth data
3.1 经典水深反演算法
3.1.1 Stumpf比值算法
Stumpf比值算法的数学表达式为
式 中，z为 水 深；
m1、m0为 模 型 的 回 归 系 数；
R(λ1)、R(λ2)分 别为蓝、绿波段的反射率；
n为蓝、绿波段的固定系数，为保证任何情况下对数值为正数且保证对数比值与深度呈线性关系而设定，本文取n=1 000。
3.1.2 Lyzenga多项式算法
Lyzenga多项式算法的数学表达式为
式中，z为水深；
a0、ai为 模型的回归系数；
N为参加反演的波段数量；
R(λi)为 波段i的反射率；
R∞(λi)为波段i对应的深水区反射率均值，可采用直方图统计方法，选择直方图中像元最小的部分作为深水区像元。
3.2 基于光谱分层的水深反演算法
本文提出基于光谱分层的浅海水深遥感反演方法，首先根据光谱对水体的穿透能力，构建光谱分层策略，提取红光波段可穿透的水深层为红光层，绿光波段可穿透但红光不可达的水深层为绿光层，蓝光波段可穿透但红光与绿光不可达的水深层为蓝光层，然后基于光谱分层构建水深反演优化模型，获取水深反演结果。
3.2.1 光谱水体穿透特性
不同光谱在水中衰减系数存在差异，在清澈水体中，蓝光、绿光、红光、近红外波段对水体穿透能力依次减弱。图2以研究区长线礁为例展示了遥感影像反射率与水深的关系，图2a中近红外波段反射率较小，且随着水深增加没有变化，说明近红外波段基本无法穿透水体；
受水体吸收的影响，图2b至图2d反射率均随着水深增加逐渐降低，其中红光波段和绿光波段分别在水深大于5 m和15 m上下时，反射率值达到最小且趋于稳定，说明红光波段和绿光波段最大穿透水深在5 m和15 m上下；
而蓝光波段在0～30 m水深区间反射率值仍未达到最小值，说明蓝光波段最大穿透水深可达30 m上下。由于水体悬浮物含量、温度、盐度、密度等因素都会影响光谱水体穿透能力，因此不同区域、不同时间的光谱水体穿透能力可能存在不同，但是蓝光、绿光、红光、近红外波段之间的水深穿透能力差异必然存在，这种不同光谱对水体穿透特性差异是本文构建光谱分层和提高反演模型精度的关键。
图2 反射率与水深散点图Fig. 2 Scatter plot of reflection values versus depth values
3.2.2 光谱分层策略
提出一种基于影像本身的无参数（无需输入任何参数）光谱分层策略。将多光谱遥感影像作为输入，利用Ostu二值化分割算法[43]，分别对近红外、红光、绿光波段进行二值化处理，提取近红外层、红光层、绿光层、蓝光层（图3）。具体步骤如下：
图3 光谱分层示例图Fig. 3 Spectral stratification example diagram
（1）近红外层Lnr提取。利用近红外波段无法穿透水体的特性，通过Ostu算法二值化近红外波段影像Inr，分割出的前景色区域为近红外层Lnr，Lnr层记录了水面高频噪声（如耀斑、船只、厚云等）与陆地区域，用于水面噪声和陆地掩膜。
（2）红光层Lr提取。利用红光最大可穿透5 m上下水深的特性，基于多光谱影像中的红色光谱影像Ir，首先，掩膜掉近红外层Lnr区域；
然后，利用Ostu算法对掩膜后的红色光谱影像Ir进行二值化处理，其前景色为水下红光信号可达区域，即红光层Lr。
（3）绿光层Lg提取。利用绿光最大可穿透15 m上下水深的特性，基于多光谱影像中的绿光波段影像Ig， 首先，掩膜掉近红外层Lnr、 红光层Lr区域；
然后，利用Ostu算法对掩膜后的绿光波段影像Ig进行二值化处理，其前景色为水下绿光信号可达但红光不可达区域，即绿光层Lg。
（4）蓝光层Lb提取。利用蓝光最大可穿透30 m上下水深的特性，基于多光谱影像中的蓝光波段影像Ib， 首先，掩膜掉近红外层Lnr、 红光层Lr和 绿光层Lg区域；
然后，将剩余区域作为蓝光层Lb。
3.2.3 分层水深反演
基于光谱分层提取的红光层Lr、 绿光层Lg、蓝光层Lb，建立基于光谱分层的Stumpf比值反演改进算法，和基于光谱分层的Lyzenga多项式反演改进算法。其中基于光谱分层的Stumpf比值反演算法数学公式如下：
式中，z为水深；
mr1、mg1、mb1分别是红光层、绿光层、蓝光层中比值模型斜率常数；
mr0、mg0、mb0分别是红光层、绿光层、蓝光层中水深0 m时，对数比值的偏移量；
n为固定常数，本文取n=1 000；
R(λr)、R(λg)、R(λb)分别为红、绿、蓝波段的反射率。此处需注意的是，不同于经典Stumpf比值算法（公式（1））中仅绿、蓝波段参与反演，在本文的改进算法中，红光层水深反演采用了红、绿波段。
基于光谱分层的Lyzenga多项式反演算法数学公式如下：
式 中，z为 水 深；
ar0和ari是 红 光 层 线 性回 归 系 数；
和是绿光层线性回归系数；
ab0和abi是蓝光层线性回归系数；
N为参加反演的波段数量；
R(λi)为 波段i的反射率；
R∞(λi)为 波段i对应的深水区反射率均值。
3.3 精度评价
为定量评价水深反演结果，采用以下3个定量指标：均方根误差（Root Mean Square Error, RMSE）、平均绝对误差（Mean Absolute Error, MAE）、平均相对误差（Mean Relative Error, MRE）。RMSE、MAE、MRE越小，水深反演效果越好。
式中，ei为反演水深值与对应的实测水深值zi之间的差值；
n为实测水深测试样本个数。
4.1 水深反演结果对比
本节将提出的基于光谱分层的反演算法与经典Stumpf比值算法和Lyzenga多项式算法对比分析。基于光谱分层的Stumpf比值算法和Lyzenga多项式算法实验设置如下：分别在红光层、绿光层、蓝光层随机选取50个水深像元作为训练样本，其余水深像元作为测试样本。经典Stumpf比值算法和Lyzenga多项式算法设置如下：上述红光层、绿光层、蓝光层3个光谱层选取的总计150个水深像元作为训练样本，同样其余水深像元作为测试样本。即所有算法共用一套训练样本和测试样本。
在长线礁实验区，对比分析经典Stumpf比值算法与基于光谱分层的Stumpf比值算法，以及经典Lyzenga多项式算法与基于光谱分层的Lyzenga多项式算法的水深反演结果（图4）。从图4来看，基于光谱分层的反演算法优于经典算法，相对于经典算法反演结果（图4a，图4c），基于光谱分层算法的反演结果在蓝色深水区斑点噪声更少，且红色浅水区域礁盘轮廓更加清晰（图4b，图4d）。图5展示了以上反演结果分别在红光层、绿光层、蓝光层和全局的误差对比，从图5整体来看，相对于经典算法反演误差图（图5a1至图5a4，图5c1至图5c4），基于光谱分层算法的反演水深与实测水深误差更聚集于0值（图5b1至图5b4，图5d1至图5d4），说明误差更小。图5中每个子图右上角误差指标可以定量地看出，基于光谱分层的反演算法明显优于经典算法，相对于经典Stumpf比值算法，基于光谱分层的Stumpf比值算法RMSE、MAE、MRE分别降低0.89 m、0.65 m、19%（图5a4，图5b4），从红光层、绿光层、蓝光层来看，误差均有降低，其中红光层MRE降低了149%（图5a1，图5b1）；
相对于经典Lyzenga多项式算法，基于光谱分层的Lyzenga多项 式 算 法RMSE、MAE、MRE分 别 降 低0.64 m、0.35 m、8%（图5c4，图5d4），从红光层、绿光层、蓝光层来看，误差同样均有降低，其中红光层MRE降低了90%（图5c1，图5d1），改进效果显著。
图4 长线礁水深反演结果图Fig. 4 Bathymetric maps of the Changxian Reef produced by inversion
图5 长线礁反演水深与实测水深值误差分布Fig. 5 Distribution of differences obtained by subtracting the measured values from the inversion water depth for the Changxian Reef
在巴克岛实验区，图6对比分析了经典Stumpf比值算法与基于光谱分层的Stumpf比值算法，以及经典Lyzenga多项式算法与基于光谱分层的Lyzenga多项式算法的水深反演结果。对比巴克岛激光测深地形图（图1c），从巴克岛右侧深水区域来看，基于光谱分层的反演算法斑点噪声较少（图6b，图6d），优于对应的经典算法反演结果（图6a，图6c）。图7展示了以上反演结果分别在红光层、绿光层、蓝光层和全局的误差对比，从每个子图右上角误差指标可以定量地看出，基于光谱分层的反演算法同样明显优于经典算法，相对于经典Stumpf比值算法，基于光谱分层的Stumpf比值算法RMSE、MAE、MRE分别降低0.41 m、0.36 m、4%（图7a4，图7b4），从红光层、绿光层、蓝光层来看，误差均有降低，其中红光层MRE降低了58%（图7a1，图7b1）；
相对于经典Lyzenga多项式算法，基于光谱分层的Lyzenga多项式算法RMSE、MAE、MRE分别降低0.45 m、0.52 m、8%（图7c4，图7d4），从红光层、绿光层、蓝光层来看，误差均有降低，其中红光层MRE降低了141%（图7c1，图7d1），改进效果显著。
图6 巴克岛水深反演结果图Fig. 6 Bathymetric maps of the Buck Island Reef produced by inversion
图7 巴克岛反演水深与实测水深值误差分布图Fig. 7 Distribution of differences obtained by subtracting the measured values from the inversion water depth for the Buck Island Reef
4.2 训练样本数量对反演精度影响
为对比分析训练样本数量N对反演算法的影响，设置不同数量的训练样本，即在长线礁实验区为15～180个，巴克岛实验区为15～600个，以15为间隔依次取值，并对于每个训练样本数量N，做100次水深反演重复实验，对100次重复实验误差值求平均作为训练样本数量为N时的反演精度，其中每次样本取值规则为：分别在红光层、绿光层、蓝光层随机等量取值，如训练样本数量为15时，则分别在红光层、绿光层、蓝光层随机取5个训练样本。
长线礁和巴克岛实验区训练样本数量对反演精度的影响如图8和图9所示，图中误差棒代表100次重复实验误差值的标准差，当标准差较大时，以数字表示，数字所在横轴位置表示对应的样本数量，数字颜色表示对应的算法。从图8和图9来看，随着训练样本数量增加，RMSE、MAE、MRE整体上呈现先降低后趋于稳定的趋势，训练样本数量N＜60时，误差呈现降低趋势，其中经典Stumpf比值算法和Lyzenga多项式算法比基于光谱分层的反演算法可以更快地实现误差稳定；
训练样本数量N≥60时，误差基本稳定且误差棒较小，该阶段中基于光谱分层的Stumpf比值算法和Lyzenga多项式算法误差低于相应的经典算法，表明训练样本数量N≥60时，可充分发挥基于光谱分层反演算法优势。
图8 长线礁实验区训练样本数量对反演精度的影响Fig. 8 Effects of the number of training samples on the accuracy of bathymetric inversion for the Changxian Reef experimental area
相对于经典Stumpf比值算法和Lyzenga多项式算法，基于光谱分层的水深反演算法拥有更高的反演精度，主要原因是基于光谱分层的水深反演算法顾及到不同光谱层的波段特性差异，红、绿、蓝波段在不同光谱层中采用了最适宜的反演系数（权重），从而取得较高的反演精度。即水深反演常用红、绿、蓝波段的测深性能存在差异，不同光谱层各波段贡献程度不一样，在红光层，因红光在水中衰减强于蓝、绿光，所以红光对水深区分度高于蓝、绿光，因此，红光层中红光对水深反演贡献高于蓝光和绿光；
在绿光层，因红光不可达，所以红光对水深反演贡献极小，因此，绿光层中蓝光和绿光对水深反演贡献最高；
在蓝光层，因红光和绿光不可达，所以红光和绿光对水深反演贡献极小，因此，蓝光层中蓝光对水深反演贡献最高。经典Stumpf比值算法和Lyzenga多项式算法全局采用一套反演系数，无法兼顾不同光谱层波段对水深反演的贡献差异，而基于光谱分层的水深反演算法根据红、绿、蓝波段极限测深能力划分光谱层，分别建立最适宜反演系数（权重），解决了上述问题，提高了水深反演性能。
从红光层MRE来看，相对于经典Stumpf比值算法和Lyzenga多项式算法，基于光谱分层的水深反演算法精度有巨大提升，MRE降低了58%～149%。主要原因如下：红光层水深较浅，极浅水域MRE反演误差的控制是当前水深反演难点，微小的水深反演偏差可引起较大的MRE，如实际水深为0.5 m的区域被反演为水深2 m，其MRE为300%。而经典Stumpf比值算法和Lyzenga多项式算法在水深较浅的红光层仍然为米级误差，造成较大的MRE。基于光谱分层的水深反演算法在水深较浅的红光层建立独立的反演系数，充分发挥红光波段优异的浅水测深性能，将红光层误差控制在分米级，大幅降低了浅水区MRE。
图8和图9显示随着训练样本数量的增加，水深反演误差整体先降低后趋于稳定，这是因为随着训练样本数量的增加，训练样本规模上可以映射出更多光谱到水深的特征规律，因此误差降低，当训练样本呈现的特征被掌握后，训练样本继续增加不会呈现新的特征规律，因此趋于稳定。其中经典Stumpf比值算法和Lyzenga多项式算法比基于光谱分层反演算法可以更快地实现误差稳定，这是因为基于光谱分层的水深反演算法对样本数量有更高的需求，即相同数量样本情况下，经典Stumpf比值算法和Lyzenga多项式算法可以使用全部样本训练出全局反演系数，而基于光谱分层的水深反演算法需要拆分为3份，分别训练对应红光层、绿光层、蓝光层的反演系数，需要更多的训练样本数据。但是随着训练样本数量增加，基于光谱分层的反演算法的反演误差迅速降低，当训练样本数量N≥60时，其误差可低于相应的经典算法，并实现基本稳定且误差棒较小，证明基于光谱分层的反演算法对训练样本数量需求并未增加太多，从水深反演应用案例来看[22-25,30-31]，N≥60的要求不会成为本文算法推广应用的障碍，训练样本数量N≥60可作为其他实验参考值。
对比基于光谱分层的Stumpf比值算法和Lyzenga多项式算法，基于光谱分层的Stumpf比值算法继承了经典Stumpf比值算法的优点，即一定程度上抑制水体不均与底质异构引起的干扰，同时也削弱太阳高度角、水面波动及卫星姿态、扫描角变化等差异对水深反演的影响[9]。基于光谱分层的Lyzenga多项式算法继承了经典Lyzenga多项式算法的优点，同样可以一定程度上抑制水体不均与底质异构问题，但相对于仅有两个波段参与运算的Stumpf比值算法，Lyzenga多项式算法可以综合多个波段的水深信息，这也是本文实验Lyzenga多项式算法精度优于Stumpf比值算法的一个重要因素（图5, 图7至图9）。总的来说，两种算法各有优势，结合实验结果与算法本身特点来看，当研究区较小且影像成像质量较佳时，即太阳高度角、水面波动及卫星姿态、扫描角变化等影响较弱时，优先推荐基于光谱分层的Lyzenga多项式算法；
当研究区较大或需要模型迁移时，有文献表明Stumpf比值算法可以取得更优的反演精度[22]。但是由于水深反演干扰因素多且海洋环境复杂，两种算法对比与适用场景值得进一步研究。
图9 巴克岛实验区训练样本数量对反演精度的影响Fig. 9 Effects of the number of training samples on the accuracy of bathymetric inversion the Buck Island Reef experimental area
本文在经典Stumpf比值算法和Lyzenga多项式算法基础上，从光谱测深性能角度，提出一种基于光谱分层的水深遥感反演优化方法。首先根据光谱对水体的穿透能力，提出一种基于多光谱本身的无参数光谱分层策略，提取红光层、绿光层、蓝光层；
然后分光谱层构建水深反演模型，建立各光谱层最适宜波段和其反演系数，解决现有水深反演算法全局一套反演系数难以适应不同光谱层的局限，最终达到提高水深反演精度的目的。以长线礁和巴克岛作为实验区开展验证分析，得出以下结论：
（1）提出的基于光谱分层的水深反演算法精度高于经典Stumpf比值算法和Lyzenga多项式算法：相对于经典算法，基于光谱分层的水深反演算法RMSE、MAE、MRE可 降 低0.41～0.89 m、0.35～0.65 m、4%～19%，证明了本文算法的有效性。
（2）在红光层，即水深较浅区域，基于光谱分层的水深反演算法MRE大幅降低：相对于经典算法，基于光谱分层的水深反演算法在红光层MRE降低58%～149%，精度提升明显。
（3）基于光谱分层的水深反演算法对样本数量与分布要求较高，但构不成算法推广应用的障碍：训练样本数量N≥60，即红光层、绿光层、蓝光层训练样本数量均大于等于20，可作为其他实验参考值。
本文从光谱测深性能角度改进了经典Stumpf比值算法和Lyzenga多项式算法并取得较好反演效果，但本文算法仍有较大改进空间。如所提光谱分层反演策略理论上可以与现有机器学习角度、地理学角度改进算法叠加使用；
如针对更多光谱的高光谱卫星数据时，通过光谱优选等方法建立更细的光谱分层策略，这些策略可进一步提升反演精度，值得未来研究。