MF多源测深数据融合方法及大洋水深模型构建

刘 洋，吴自银，赵荻能，周洁琼，尚继宏，王明伟,朱 超，鲁号号

1. 自然资源部第二海洋研究所，浙江 杭州 310012； 2. 国家海洋局海底科学重点实验室，浙江 杭州 310012； 3. 上海交通大学海洋研究院，上海 200030

精确可靠的水深数据几乎是从事一切海洋活动的基础[1]，高分辨率的数字水深模型(digital bathymetric model，DBM)不仅可用于海图的构建，还可用于地貌特征与成因[2-3]、海底构造、海洋资源、水动力流体模型、生物栖息地与生态系统等多个方面的研究[4-14]。对海洋测深的系统研究持续了一个多世纪[15-19]，从早期极低效率的锤线测深发展到当前高精度、高分辨率、高效率的多波束测深技术[20-22]。然而，由于全球海洋极其广阔，基于船载模式的多波束测深成本高昂[23]，导致目前仍有超过80%的全球海域未被多波束探测[24]。因此，在构建深海大范围DBM时仍需使用历史测深资料[25]。当前大洋测深数据主要针对不同目的独立采集，导致来源、格式和分辨率等方面存在诸多差异[26-28]。全球大洋测深数据的多源复杂性，不利于高精度DBM的构建，更不利于日益迫切认识海洋、经略海洋的国家战略需求。

DBM的构建也称数据的网格化。网格化是二维采样、插值和二维外推的结合[29]。如何基于插值算法将多源水深数据融合进而构建DBM一直是学者们的研究热点。根据目前国内外对多源水深数据融合技术的研究，融合方法主要是基于克里金、双三次样条、张力样条等插值算法进行综合应用，如文献[30]将多波束、单波束、历史海图等水深数据进行融合，编制了北冰洋的国际海底地形图(international bathymetric chart of the arctic ocean，IBCAO)；文献[31]根据多种水深数据汇编了南大洋的国际海底地形图(international bathymetric chart of the southern ocean，IBCSO)；文献[32]结合不同时期的多波束、单波束和海岸线等数据，构建了南极洲东部George V大陆边缘的DBM，得到该区较详尽的海底地貌信息。上述方法虽实现了多源水深数据的融合，但融合后的数据空间分辨率低(>500 m)，且难以保证数据稀疏区和空白区的插值精度。因此，在保证数据稀疏区插值精度的同时，既要提高数据空白区的精度，又要保留高分辨率水深源数据的细节特征，是值得深入研究的科学技术问题。

为解决上述问题，本文提出合并-融合(merge-fusion，MF)法。该方法将多波束、单波束、电子海图数据与通用全球海洋地形(general bathymetric chart of the oceans，GEBCO)数据融合，在保留高分辨率地形细节特征的同时，合理填补了数据空白区，形成“挑战者深渊”分辨率100 m的高精度DBM。结果表明，该方法融合的DBM能更好地反映精细的地形特征信息，为区域性或全球性DBM的构建以及多源异构水深数据的管理利用提供了实用化技术方案。

1 MF方法基本原理与步骤

MF(merge-fusion)法的核心思想为“合并-融合”。如图1所示，该方法主要包括6个步骤：

(1) 原始数据处理与入库。根据不同水深数据的来源、类型，对其进行精细处理。以多波束数据为例，对于系统误差，通过潮位改正、声速剖面改正、换能器安装偏差校准等方法进行消除[33-39]；对于异常值，通过人工或者自动滤波算法识别并剔除[20]；而随机误差，只能通过对测深数据的精度评估评价其误差大小[40]。建立多源水深数据库，对清理后的水深数据进行统一管理。

(2) 多源数据精度评估。对多源水深数据进行精度评估是保证成果网格数据质量的重要手段[40]。在同一水域，以多波束条带数据的中央波束测深结果为基准，使用平均差值式[30](1)对多源数据进行精度评估

(1)

式中，δ为平均水深差；q为测深数据个数；DMi为多波束水深数据；DSi为多源水深数据，水深值取负值。

根据《GB 12327—1998 海道测量规范》，将水深偏差与水深的2%比较，若偏差小于水深的2%，符合规范要求不做处理；否则根据偏差值以多波束条带数据的中央波束测深结果为基准，对其他类型数据进行最小二乘法平差处理[28,41]。通过式(2)来表示误差效应

F(x,y)=a0+a1x+a2y+a3x2+a4xy+a5y2

(2)

式中，(x,y)为测点坐标；a0,a1,a2,a3,a4,a5为待定系数。

设水深观测值Z的改正数为V，以多波束条带数据的中央波束测深结果为基准，建立如下误差方程

-(VSi-VMi)=(a0+a1x+a2y+a3x2+

a4xy+a5y2)Si-(a0+a1x+

a2y+a3x2+a4xy+a5y2)Mi-

(ZSi-ZMi)

(3)

式中，VSi、VMi分别为多源数据和多波束数据在点i位置(x,y)处的改正数；ZSi，ZMi分别为相应位置处的多源数据和多波束数据水深值。

图1 MF方法原理及流程Fig.1 Principle and workflow of the MF method

根据∑(VSi-VMi)2=min解算误差方程(3)得到误差改正系数，最后依据式(2)对所有测深数据进行改正，即获得平差后水深值。

(3) 多源数据网格化。为避免网格内潜在异常值的影响[42]，对多源水深数据库导出的离散水深点进行中值滤波。中值滤波首先确定奇数像元网格G，网格内各像元按数值大小排序，并将网格节点内的中值水深d(x,y)作为原数值s(x,y)[43]

d(x,y)=median{s(x-u,y-v),(u,v)∈G}

(4)

为保证数据稀疏区域插值的准确性[44]，对滤波后的多源水深数据进行网格化[45-47]，形成低分辨率网格。地球科学领域最常见的网格化方法之一是双三次样条插值算法[48]，该算法可通过已知点精确插值，得到光滑的插值曲面，但在地形突变处或数据稀疏区易产生异常值，掩盖真实地形。本文采用目前国际上普遍使用的连续曲率张力样条插值算法[48-49]对多源数据进行网格化，该算法由于引入张力因子，可消除最小曲率拟合面中存在的较大波动和无关变形点[43]，网格化后的结果更加平滑、准确。该算法基于最小曲率拟合面原则

(5)

式中，(xi,yi,zi)为已知点坐标；N为已知点个数；fi为点(xi,yi)处格林函数线性组合解中的系数；δ(x-xi,y-yi)为给定响应函数；2为拉普拉斯算子。

边界条件为

(6)

在此基础上，引入张力因子T，得到

(1-TI)2(2z)-TI2z=

(7)

式中，TI为拟合面内部张力因子，在[0,1)区间内取值。

边界条件为

(8)

式中，TB为拟合面边界处张力因子，在[0,1)区间内取值。

当张力因子T=0时，式(7)简化为式(5)，此时无张力，得到最小曲率平滑解；T=1时，张力最大，得到调和样条解。对于大洋水深数据，张力因子一般选取0.32～0.4之间[30]。随后对低分辨率网格滤波，将滤波后的低分辨率网格重采样至目标分辨率，得到“基础水深网格”。

(4) 多分辨率网格合并。为合并多分辨率网格，且保留高精度测深区域地形细节特征，使用高分辨率网格水深值替换低分辨率网格水深值。首先，将预处理步骤中输出的高分辨率网格与基础水深网格叠加分析，得到二者在同一位置处z的差值；对z的差值进行判断，若差值非零则用多波束数据网格点内相应值进行替换，得到“合并水深网格”。

为消除网格合并过程出现的边缘效应[50]，在高分辨率网格数据一侧定义过渡带，过渡带宽度至少为3列网格点，具体范围根据接边情况来确定。其原理基于双曲线加权函数

(9)

式中，L代表到下一高分辨率数据约束单元或过渡带外约束单元的距离；w表示权重[31]。

由式(9)可看出，距离高分辨率数据网格单元越近，受深度影响越小，由此弱化边缘效应发生的概率。

(5) 数据空白区融合。在数据空白区域范围内，将GEBCO网格与合并水深网格“作差”，即GEBCO网格内每个数据点减去相同位置处插值数据点的值，得到包含经度、纬度和水深差异的差值文件。随后网格化中值滤波后的差值文件，得到差值网格。将差值网格叠加整合到合并水深网格，得到融合后的DBM。

(6) 网格可视化及清理。经上述步骤得到的DBM难免存在各种“毛刺”等异常点，使用三维地形浏览法找出这些问题。为突出数据噪声，可使用光照效果加以突出，最后通过曲面拟合算法[51]对网格问题区域进行编辑和误差处理，以获得最佳的数字水深模型。曲面拟合算法的构建原理：根据波束点拟合海底趋势面，计算实测水深数据与趋势面间的深度差值，结合误差处理理论剔除异常点[52]。曲面拟合法中拟合函数的一般形式

(x,y,z)∈Q(xc,yc,zc)}

(10)

式中，Ψ(x,y)为曲面拟合函数；(x,y,z)为波束点空间坐标；hlm为多项式系数；k为多项式总阶数；(xc,yc,zc)为待检测点c；Q(xc,yc,zc)为以c为中心的曲面拟合函数的局部曲面拟合范围。

2 数据来源

位于马里亚纳海沟南端的“挑战者深渊”是全球海洋最深处，被誉为“世界第四极”，一直以来都是全球海洋科学研究的热点区域[53-56]。近几年来，我国在该区域组织了数次综合海洋调查。高精度的海沟区DBM是开展载人深潜、海底原位观测和板块俯冲数值模拟[57]的基础，但该区现有数据分辨率低、精度差，难以满足海洋科学研究和工程实施要求。同时，该海区数据来源多样，且存在诸多数据空白区，适合数据融合技术的方法研究。

本文数据主要来自于公开的海洋数据，包括多波束数据、单波束数据、电子海图数据以及GEBCO数据。

2.1 多波束数据

构建此次DBM的核心数据主要是多波束数据，采集自2003年以后，且精度较高。在该区域共收集25个航次的多波束数据，其中，15个调查航次数据源于美国大气与海洋局(NOAA)下属的国家地质与地球物理数据中心(NGDC)，另外10个多波束测深航次数据则从日本海洋-地球科学和技术研究中心(JAMSTEC)获得。

2.2 单波束数据

除多波束数据外，还从NGDC收集了2000年以来的大量单波束测深点数据(图2)。由于收集的单波束数据未提供测深仪及定位系统等相关详细信息，很难对所有数据进行精度评估。且单波束测深点空间分布非常不均匀，只在使用多波束探测区域分布较为密集，例如美国关岛附近。

2.3 电子海图数据

从美国国家海洋和大气管理局(NOAA)下的海岸调查部门(Office of Coast Survey)共收集该区域2005年来的4幅电子海图数据。除电子海图数据外无相关元数据信息，因此无法全面评估其数据精度。

G区为多源数据精度评定所选区图2 “挑战者深渊”位置和数据源Fig.2 Location and data sources of the “Challenger Deep”

2.4 GEBCO数据

虽然已有大量水深数据，但仍存在范围达上千平方公里的数据空白区，本文使用GEBCO_2014[58]水深数据集对数据空白区进行填补。GEBCO数据是当前国际主流水深数据集之一，而GEBCO_2014为目前全球测深数据集的最新版本，主要源于实测数据和利用卫星高度计反演的预测水深，网格分辨率为30″。

3 MF方法的实现

根据上述方法对实验区多源水深数据进行融合处理，构建马里亚纳海沟“挑战者深渊”的DBM。首先将多源水深数据处理为离散数据格式，转换水平坐标，投影方式和深度基准面分别统一至UTM投影和平均海平面。随后进行系统误差处理，并使用CUBE算法清理多源水深数据的异常点[59-60]，后重新导出为离散格式，形成“挑战者深渊”的多源水深数据库。

为对多源水深数据进行精度评定，比较了美国关岛西北部区域(图2区域G)内多源水深数据的深度差异。该区域被多波束数据全覆盖，精度较高，适合用于精度评估。离散的电子海图和单波束水深点分别与分辨率100 m的多波束条带数据的中央波束测深结果进行比较(图3)。两者分别与多波束进行一一组合，组合的准则是二者间距小于50 m[30]。选取115对电子海图水深对比数据、207对单波束水深对比数据，依据式(1)得到对比结果，分别见表1和图3。

蓝色点代表电子海图数据，红色点代表单波束数据，直方图表示深度差值绝对值分布图3 电子海图、单波束与多波束数据水深值对比Fig.3 Comparison of water depth of electronic charts, single-beam and multi-beam bathymetric data

Tab.1 Comparison of multi-source dataaccuracyassessmentm

对象平均差值平均水深±2%·水深电子海图vs多波束-18.9-1039±20.7单波束vs多波束-23.2-1278±25.5GEBCO_2014vs多波束-10.7-1185±23.7

根据《GB 12327—1998 海道测量规范》，对于电子海图数据，平均水深1039 m时，平均水深差值优于水深值的2%；对于单波束数据，平均水深1278 m时，平均水深差值同样优于水深的2%。经过评定，单波束数据、电子海图数据精度在限差容许范围内，可作为本次数据融合实验的数据集。

以区域G内多波束条带数据的中央波束测深结果为基准，对GEBCO_2014数据进行精度评估。在分辨率100 m基础上，得到二者水深差值的空间位置分布(图4)。应用式(1)从29605对GEBCO对比数据中，得到对比结果，分别见表1和饼状图(图4)。平均水深1185 m时，平均深度差优于水深的2%。GEBCO_2014数据精度同样在限差容许范围内，可作为本次数据空白区融合实验的数据集。

图4 GEBCO_2014数据与多波束数据水深值对比Fig.4 Comparison of water depth of GEBCO_2014 data and multi-beam bathymetric data

对多源水深数据(图5(a))进行网格大小1000 m的中值滤波。随后网格化滤波后的水深数据，形成分辨率1000 m的低分辨率水深网格(图5(b))。对低分辨率网格进行平滑和滤波，将滤波后的低分辨率网格重采样至100 m作为“基础水深网格”(图5(c))。

通过对多波束数据进行处理，例如声速剖面改正、异常值剔除等测深条带拼接步骤[33,43,61-63]得到分辨率100 m的高分辨水深网格(图5(d))。然后将高分辨率水深网格与基础水深网格叠加分析，得到二者在同一位置处z的差值。若某点处z的差值非零，则用高分辨率水深网格内相应水深值替换，得到“合并水深网格”(图5(e))。为消除多分辨率网格间的边缘效应，在靠近高分辨率网格数据边界定义5 km的过渡带(图6)，极大地降低了边缘效应产生的概率。

图5 多源数据融合过程示例Fig.5 An example of multi-source data fusion

图6 过渡带定义前后对比Fig.6 Comparison of the DBM

结合GEBCO_2014数据，对数据空白区域进行融合。首先将GEBCO_2014网格与合并水深格网叠加合并，得到数据空白范围内同一位置处的水深差值文件，差值文件包括经度、纬度和深度差值信息。随后对差值文件进行网格大小1000 m的中值滤波。将滤波后的差值文件网格化为分辨率100 m的差值网格。最后，将差值网格与合并水深网格叠加整合，得到融合后的DBM(图5(f))。

可视化浏览合并-融合后的水深网格，识别并高亮异常点，对网格问题区域进行编辑和误差处理，最终得到满意的DBM。

4 分析与讨论

4.1 数据拼接对比

图7为使用MF法得到的马里亚纳海沟“挑战者深渊”DBM，网格分辨率100 m。经统计，本次所选区域多波束、单波束和电子海图数据的覆盖率分别为42.2%、0.5%、0.2%，使用GEBCO数据对剩余57%的数据空白区进行填补。通过网格可看出，MF方法不仅保留了高分辨率数据区详细的海底特征，还为数据稀疏区和空白区提供了平滑的地形表面。

在马里亚纳海沟“挑战者深渊”地区，除本文构建的“挑战者深渊”DBM之外，只有全球数据集涵盖整个大洋。本文将“挑战者深渊”DBM与GEBCO_2014DBM进行比较。结果表明，在数据接边处，MF法取得较为理想的效果，接边处总体限差符合《GB 12327—1998 海道测量规范》，即水深大于100 m时，深度测量的极限误差应优于水深的2%。在实验区选取接边差异较大区域M(图7)，针对MF法和传统方法得到剖面a-a′(图8(a))和剖面b-b′(图8(b))。从剖面图可看出，两种方法形成差异主要集中在距剖面线起点2～4 km范围内区域。以剖面线起点为原点，剖面方向为正方向，选取(2.6 km，2.8 km)范围内剖面A和B(图8)，在水深约3600 m处，水深限差50 m时，对于传统方法，平均坡度为72.32°(图8(b))；而MF法，平均坡度为37.63°，与前者相比整体过渡较为平缓(图8(a))。

4.2 数据细节对比

高分辨率水深数据结合越多，网格分辨率就会越高，这使得本次融合的分辨率100 m的“挑战者深渊”DBM成为该区域目前精度最高的数据集；同时，该数据集对“挑战者深渊”地区的地形细节特征刻画有着显著改善。例如，在网格单元100 m的“挑战者深渊”DBM中，可清晰地看到海底峡谷等地形特征(图9(d))，以及形态、轮廓更明显的海山(图9(c))，而在GEBCO_2014数据集中几乎看不到此类地形特征。在数据稀疏区，利用清理过的单波束或电子海图数据形成的低分辨率网格避免了异常点的出现(图9(a)和9(b))。但仍存在数据密度和质量较差地区，这些区域包含潜在异常点、噪声点，需进一步调查才能识别和处理。

5 结 论

(1) 提出并实现了一种适用于多源水深数据融合的MF方法，并给出详细的融合技术方法流程。其基本原理是基于多源水深数据形成两种分辨率网格，然后将多分辨率网格合并，在此基础上使用GEBCO数据进行融合，叠加、合并得到最终的DBM。

M为数据拼接对比区；a、b、c、d 为数据细节对比区图7 “挑战者深渊”DBMFig.7 The DBM of the “Challenger Deep”

图8 数据接边处比较Fig.8 Comparison of data transition area

a1、b1、c1、d1为MF法所得DBM，a2、b2、c2、d2为传统方法所得DBM，A代表异常点，C代表海底峡谷图9 “挑战者深渊”DBM与GEBCO_2014DBM细节对比Fig.9 Comparison of data details between “Challenger Deep” DBM and GEBCO_2014 DBM

(2) 基于马里亚纳海沟“挑战者深渊”地区公开的多波束、单波束、电子海图等水深数据，采用MF多源数据融合方法，构建了“挑战者深渊”的高精度DBM，验证了本文方法的可靠性和实用价值。

(3) 针对同一区域，利用本文方法得到的DBM与目前国际主流的GEBCO_2014数据集进行对比，表明本文方法更能够体现高分辨率水深数据和数据空白区的细节信息、防止边缘效应和异常点的产生，从而获得高精度的海底DBM。