本文内容来源于《测绘学报》2023年第12期（审图号GS京(2023)2441号）

基于最小二乘配置的陌生海域海底地形反演方法

范雕1, 李姗姗1, 冯进凯1, 黄炎2, 范昊鹏1, 张金辉1, 李新星1

1. 信息工程大学, 河南 郑州 450001;

2. 军事科学院国防科技创新研究院, 北京 100071

基金项目：国家自然科学基金(42204009;42174007;42174008)

摘要：针对目前在无船测水深数据支撑的陌生海域海底地形模型构建困难的现状, 提出了可将船测水深数据和重力数据已知海区协方差函数统计结果迁移、推广并应用于与已知数据海区重力异常相似的船测水深数据空白的陌生海区, 运用最小二乘配置方法实现陌生海域海底地形模型构建。首先在西太平洋某1°×1°数据已知海区, 依据最小二乘配置方法构建了相应的海底地形模型(BAT_LSC_1)。模型评估结果表明, BAT_LSC_1模型与导纳函数构建的海底地形模型检核精度相当, 检核相对精度优于4%。海区重力异常作灰度化处理后, 应用影像粗匹配和精匹配技术搜寻出与已知海区重力异常特征相似区域, 迁移已知海区协方差统计结果应用于陌生海区, 运用最小二乘配置方法构建陌生海区海底地形模型(BAT_LSC_2)。试验检核结果表明, BAT_LSC_2整体表现优于ETOPO1模型和应用导纳函数构建的海底地形模型, 验证了本文方法思路的可行性和适用性。

关键词：海底地形 卫星测高 重力异常 最小二乘配置 协方差函数

引文格式：范雕, 李姗姗, 冯进凯, 等. 基于最小二乘配置的陌生海域海底地形反演方法[J]. 测绘学报，2023，52(12)：2039-2053. DOI: 10.11947/j.AGCS.2023.20220644

FAN Diao, LI Shanshan, FENG Jinkai, et al. Applying least square collocation method to predict seafloor topography in the unknown sea area[J]. Acta Geodaetica et Cartographica Sinica, 2023, 52(12): 2039-2053. DOI: 10.11947/j.AGCS.2023.20220644

阅读全文：引 言

海洋是人类生存和发展的重要资源宝库，其面积占地球表面积71%左右[1]。文献[2]研究表明目前依据船载声学测量手段仅有15%左右的海域面积完成了空间分辨率约800 m的海底地形绘制，人类对于海洋的认识仍处于较低水平和层次[3]。文献[4]认为依托船舶为载体实施海底地形测量，若以空间网格大小为1′作为考虑对象，全球大约还有90%的海域属于未被勘探状态，而且船舶测线间存在巨大的空间间隔，测量结果分布很不均匀，人类对海底地形的描绘详细程度远不及对月球等其他行星的了解深入[5]。

然而，海底地形地貌数据作为开发和利用海洋的基础性地理空间信息数据，广泛影响其他科学研究和应用实践。一方面，海底地形变化直接影响海洋环流和混合，从而间接影响地球气候、生物多样性和粮食产量。另一方面，海底地形和海底构造的精细化表达有利于研究板块构造规律，有利于揭示海盆和海底地质运动历史，有利于预判发生地震、火山爆发和海啸等区域位置。因此，过去100多年来国内外学者持续不断努力研究绘制全球范围海底地形图[6-7]，囿于技术装备发展水平制约，前期研究进展较为缓慢。直到19世纪90年代中期，由于GEOSAT和ERS-1等执行大地测量任务的测高卫星数据不断积累和应用，全球海洋区域重力场精度和分辨率不断提升，依托海洋重力信息反演海底地形研究不断推进[8-10]。至此，运用空天探测技术获取全球范围的海底地形信息方才取得重要突破。文献[11]利用卫星测高重力异常数据，依据回归方法(又称S&S方法)建立了范围涵盖南北纬72°之间，分辨率约15 km的海深格网模型。随后数十年时间，由于模型构建数据量不断丰富，数据精度不断提升，学者们在更高精度、更高分辨率及更广覆盖范围的海底地形建模方面进行了积极探索，比如：①丰富的船测海深数据可用于校准和修正预测海深结果[12]，以使构建的海深模型网格密集程度更高；②波形重跟踪技术的广泛使用，在算法层面不断提高测高结果精度，从而不断提升海洋重力场质量[13]。目前利用重力数据反演海底地形算法大致可分为统计法和解析法两大类：统计法包括回归分析方法[14]、模拟退火方法[15]、最小二乘配置[16]及神经网络[17]等；解析方法包括重力地质方法[18-19]、导纳函数方法[20]及顾及地形高次项影响的解析反演方法[21-22]等。然而，重力地质方法、回归分析方法、模拟退火方法在试验海区需要外部船测水深数据支撑，在无船测水深数据支撑的陌生海区该方法将不再适用，进而极大削弱算法应用范围。基于Parker公式的解析反演方法不需要外部船测水深数据为补偿，文献[23]利用此方法开展了无船测水深数据(利用先验海深模型代替)的海底地形建模，但是获取该方法涉及的诸多地球物理参数并非易事，若全球海域均采用参数经验值代替，毋庸置疑会造成最终反演精度损失，进而导致反演结果可信度降低。

基于此，为扩展无船测水深数据支持环境下陌生海域海底地形恢复方法，本文提出以卫星测高重力数据和先验海底地形模型为输入，运用最小二乘配置方法(统计法)，解决无船测数据支撑环境下高效构建陌生海区海底地形模型难题的新思路。基本策略是：首先，在船测水深数据和重力数据已知海区统计解算出反映海区地形的相关协方差模型；然后，采用最小二乘配置方法，以海区重力数据为输入，构建相应的海底地形模型并检验算法适用性；最后，将已知海区数据统计结果(不同类型参量的协方差函数)迁移、推广并应用于与已知海区相似的船测数据陌生海区，实现仅依靠重力数据的海底地形模型构建。同时将本文反演结果与导纳函数反演海底地形结果、DTU18、ETOPO1模型进行比对分析，简要评估模型精度和效能。

1 原理与方法

1.1 最小二乘配置方法

利用解析方法(例如导纳函数)反演海底地形时，算法模型涉及的地球参数求解并非易事。文献[24]认为海底地形和重力信息为两类随机变量，从而依据重力异常数据，采用统计方法恢复海底地形模型为

(1)

式中，b为海底地形；g为海面重力异常信息；μb和μg表示海底地形和海面重力异常的数学期望(可用相应的平均值代替)；Cbg为海底地形和重力异常的互协方差矩阵；Cgg为海面重力异常的自协方差。式(1)化简为

(2)

式中，b为b－μb；g为g－μg。式(2)即为最小二乘配置方法表达式[25]。

由式(2)可知，使用海面重力异常数据，运用最小二乘配置方法恢复海底地形信息关键在于获取地形-重力异常互协方差函数和重力异常信息间自协方差函数。其中重力信息的自协方差函数求解在物理大地测量重力拟合推估研究较多[26]，目前已有成熟的解决手段。然而海底地形和重力异常分属于两种不同类型的参数，因此如何精确求解地形-重力异常互协方差函数成了利用式(2)恢复海底地形的主要难题。

1.2 地形-重力协方差估计

将地形和重力异常作为一维随机变量参考时，地形-重力异常互协方差的无偏估计通过遍历数据得到[27]

(3)

式中，g(x)和b (x)表示经过归中处理后的重力异常和地形；T为数据周期；ξ为地形-重力异常间距离。可见，重力-地形互协方差函数为距离ξ的函数。需要注意的是，一般易证明，它们之间的关系为[27]

(4)

式(3)离散化的结果为

(5)

式中，l表示重力异常和地形间距离，通常取0，1，…，N-1。易知，式(5)重力-地形协方差估计是无偏估计

(6)

式中，E(·)表示期望算子。同理，二维重力-地形互协方差函数无偏估计表示为

(7)

借鉴重力异常拟合推估自协方差函数求解策略[26, 28]，对于方形格网研究区域(N行和M列)重力-地形的互协方差可采用式(8)统计纵向和横向互协方差

(8)

式中，q和qj分别表示横向距离和纵向距离。最终将所有具有距离l的纵向协方差和横向协方差取平均值作为协方差Cbg(l)的统计数值

(9)

协方差统计数值结果反映了海区相关数据变化的内在规律，通常情况下需构建强相关函数模拟海区数据协方差统计结果，进而依托强相关函数拟合结果构造相应的(互)协方差矩阵，从而据最小二乘配置算法模型恢复海底地形信息。强相关函数如Inverse协方差函数、Gauss协方差函数、Hirvonen协方差函数、似高斯协方差函数、周期指数函数、Tscherning/Rapp协方差函数和Mortiz协方差函数等[29]。本文选择较为简单的似高斯协方差函数作为海区协方差统计备选拟合函数，表达式为

(10)

式中，C0是两点距离为零对应的协方差函数值，表示协方差函数的变化幅度；l为两点之间距离；a为常数。式(10)通过恒等变换后为

(11)

从而依据最小二乘方法可求解似高斯协方差函数待求参量a。借鉴以上分析结果，同理可求得重力-地形协方差 Cgb和重力异常的自协方差 Cgg，从而可获取以海面重力信息为输入，采用最小二乘配置方法(式(2))构建相应波段的海底地形模型。

2 数值分析试验

2.1 数据准备与前期处理

选择西太平洋某1°×1°(155.2°E—156.2°E，22.5°E—23.5°E)海区作为研究试算区域，海区输入数据包括卫星测高重力异常数据和稀疏船测数据，数据来源和预处理如下。

(1) 海区重力异常来源于SIO(Scripps Institution of Oceanography)于2022年8月3日发布的S&S V32.1重力异常模型，该版本重力异常模型增加了Altika、Cryosat LRM、Cryosat SAR和Sentinel-3A/B等数据。研究海区重力异常空间分布如图 1所示，相应的统计结果见表 1。重力异常较大值集中在海区中部区域，推测该区域可能存在明显的海山地形。

图 1 研究海区重力异常Fig. 1 Gravity anomaly in study sea area

图选项

表 1 研究海域重力异常(S&S V28.1)数据统计Tab. 1 Gravity data statistics of research sea areamGal

表选项

(2) 稀疏船测海深数据来源于NGDC(The National Geophysical Data Center)发布的单波束和多波束船测海深数据。共收集到1061个实测船测点，数据频数分布直方图如图 2(a)所示。依图 2(a)可知，研究海区水深测量结果主要集中于5000 m以上。采用一定的处理手段对原始水深测量数据进行质量控制，本文处理策略是依据2019年12月发布的海洋空间格网大小为15″的SRTM15+V2.0作为参考模型[30]，采用双线性内插技术将SRTM15+V2.0模型海深值内插到离散船测点处，并与船测点海深值作差得到海深残差；借鉴3σ误差处理方法，认为海深残差绝对值大于3倍海深残差标准差的船测点是粗差点，然后作剔除处理。最终利用3σ方法处理原始船测海深数据后，共剔除34个可疑水深点，剩余1027个船测水深点。原始水深、处理后水深及水深残差统计结果见表 2，表中括号内数值代表参与统计的水深点个数。由表 2可知，以SRTM15+ V2.0作为对比模型，海深残差标准差为44.59 m，从而粗差阈值为133.77 m。原始船测水深数据剔除粗差后结果的相应直方图如图 2(b)所示，船测数据空间分布情况如图 3所示(图中背景为S&S V19.1海深模型)。

图 2 船测数据分布直方图Fig. 2 Distribution histogram of shipborne bathymetry data

图选项

表 2 海深数据统计结果Tab. 2 Statistical results of bathymetric datam

表选项

图 3 研究海区粗差剔除后船测点分布Fig. 3 Distribution of survey points after gross error elimination

图选项

利用卫星测高重力数据通常仅可恢复有限波段海底地形信息，剩余频段地形信息常使用实际船测水深数据进行补偿；另外为评估最终构建的海深模型精度，通常选择部分船测点作为外部检核参考[21]。从而本文基于1027个粗差处理后的船测点，随机均匀选择847个船测点(约占船测点总数的82.47%)作为控制点以补偿非反演频段地形信息；剩余180个船测点(约占船测点总数的17.53%)作为评估海深模型的外部检核点。

另外，S&S V19.1海深模型中奇数海深值由实际海深约束生成，偶数海深值为重力反演预测水深值，从而可认为这些奇数海深值是实际船测海深值。研究海区S&S V19.1海深模型中奇数(实际船测数据约束)和偶数(重力预测海深)海深值空间分布情况如图 4(a)所示，其中“·”代表重力预测海深的位置，“ ”代表实际水深约束生成的海深值位置。S&S V19.1模型在试验海区共有616个船测实际值，相应的统计结果见表 2，海深平均值和标准差分别为4 861.00 m和816.09 m；船测数据频数分布直方图如图 4(b)所示。

图 4 S&S V19.1船测水深数据Fig. 4 S&S V19.1 shipborne bathymetry data

图选项

对比S&S V19.1模型中616个船测点与笔者收集到的1027个船测点空间分布如图 5(a)所示，图中“·”号表示笔者收集的船测点位置，“”号代表S&S V19.1中船测点位置。图 5(a)船测点空间分布比较就总体而言，笔者收集的船测数据与S&S V19.1在大部分区域重合；部分海区S&S V19.1相比于本文增加了多条测线的数据，如研究海区北部，S&S V19.1中包含了笔者不掌握的船测水深数据。需要说明的是，因S&S V19.1空间网格大小为1′，从而基于S&S V19.1提取的船测点间距也是1′；而笔者从NGDC下载的声学水深测量结果为不规则离散分布，故认为S&S V19.1使用了本文以外的其他水深测量数据，数据量更加丰富。从而为尽可能充分利用可用船测数据，将S&S V19.1中包含的船测点(616个)和本文收集的船测控制点(847个)合并以补偿非反演频段海底地形信息(共1463个)。最终，试验海区控制点和检核点空间分布如图 5(b)，其中三角形表示检核点位置，白色六边形表示控制点的空间位置。

图 5 船测水深点分布Fig. 5 Distribution of shipborne bathymetry points

图选项

2.2 海底地形模型构建

海底地形主要影响海面重力中长波段信息，因而利用卫星测高重力数据也主要恢复海底地形中长波段信息。依照文献[31]和笔者前期研究相关结论，本文选择15 km和160 km作为截止波段，即依托卫星测高重力异常数据主要恢复15~160 km波段海底地形信息，15 km以内和160 km以外的海底地形信息依托船测海深进行补偿。首先，利用GMT(generic mapping tools)中格网化操作模块(surface)将试验海区离散船测控制点格网化为网格大小为1′的海深格网[32]；然后，利用文献[21]基于挠曲补偿和信噪比设计的高通、低通和带通滤波器对格网海深和S&S V32.1重力异常实施滤波操作，分别获得长波(波长大于160 km)和短波(波长小于15 km)海底地形，反演波段(15~160 km)重力异常，相应的数据空间分布和统计结果如图 6和表 3所示。需要说明的是，依据海底地形解析反演算法[33]可知，重力异常经带通滤波并向下延拓后与海底地形表现出较强的相关关系，因而为充分利用重力场元与海底地形间相关关系进行统计分析，图 6(c)表示的反演波段重力异常为经带通滤波并延拓到平均海深面的重力异常结果。

图 6 不同波段数据Fig. 6 Different band data

图选项

表 3 数据统计结果Tab. 3 Statistical results of data

表选项

对比图 6(a)和图 6(b)可发现，长波海底地形表现一种大尺度的低频特征，而短波海底地形数据则呈现出高频变化特征，从而可利用重力梯度和重力异常对不同频段地形敏感程度进行联合反演[34]。依据式(2)恢复海底地形信息需解决的关键技术问题是地形-重力协方差矩阵和重力协方差矩阵的获取，其中重力异常协方差矩阵为对称矩阵，地形-重力协方差矩阵是非对称矩阵。重力异常协方差矩阵求解技术目前已十分成熟[35]，而地形-重力协方差矩阵研究成果相对较少。分析地形-重力协方差矩阵元素结构，结合互协方差性质易知，地形-重力协方差矩阵中上三角矩阵元素值可将对应地形和重力位置关系施以地形-重力协方差函数(如本文选择的似高斯协方差函数)得到；地形-重力协方差矩阵下三角矩阵元素值则是对应重力和地形位置关系施以重力-地形协方差函数获取。

基于以上分析，首先以反演波段重力异常格网数据为输入，统计分析试验海区重力异常分布规律，为使统计结果充分反映统计数据间相互关系，参与统计样本应尽可能丰富，从而本文实施数据统计分析时，数据间横向距离和纵向距离最大设置为40′。考虑到协方差函数非负性质，去掉协方差统计为负值结果，最终得到的重力异常协方差统计结果如图 7(a)所示，距离为零时，协方差统计结果为2 721.248 2 mGal2。依据重力异常协方差统计结果，根据式(11)，采用最小二乘方法拟合参数结果如图 7(b)所示，图中方点为协方差待拟合数值分布，直线为拟合结果，相应的参数计算结果为0.207 8。采用相同的统计策略，得到的重力-地形互协方差和地形-重力互协方差结果如图 8和图 9所示。其中距离为零时，重力-地形互协方差和地形-重力互协方差统计值大小均为33.776 7 mGal·km，其中重力-地形互协方差似高斯函数参数拟合值为0.210 8，地形-重力互协方差似高斯函数参数拟合值为0.192 1。最终，基于似高斯协方差函数的(互)协方差函数表达式为

图 7 重力异常协方差结果Fig. 7 Gravity anomaly covariance results

图选项

图 8 重力-地形协方差结果Fig. 8 Cross covariance results between gravity and terrain

图选项

图 9 地形-重力协方差结果Fig. 9 Cross covariance results between terrain and gravity

图选项

(12)

研究海区格网点间欧氏距离(距离单位为′)解算结果如图 10所示，对图 10距离矩阵施以式(12)重力异常似高斯协方差函数解算，可得相应的重力异常协方差矩阵(如图 11(a))。依前面分析知，地形-重力协方差矩阵为非对称矩阵，从而将距离矩阵所有元素直接施以地形-重力协方差函数计算结果并非地形-重力协方差矩阵。对于地形-重力协方差矩阵求解，分为两个步骤：①对图 10距离矩阵的上三角所有元素施以地形-重力似高斯协方差函数，距离矩阵的下三角矩阵所有元素施以重力-地形似高斯协方差函数；②将各自计算结果合并可得地形-重力协方差矩阵，结果如图 11(b)所示。基于获取的试验海区地形-重力协方差矩阵和重力异常协方差矩阵，以反演波段重力异常为输入，采用最小二乘配置方法可恢复相应波段的海底地形，然后叠加依据稀疏船测数据获取的非反演波段海底地形可获得试验区域海深模型。试验海区海深模型结果如图 12所示，为描述方便，称为BAT_LSC_1模型。

图 10 距离矩阵Fig. 10 Distance matrix

图选项

图 11 协方差矩阵Fig. 11 Covariance matrix

图选项

图 12 BAT_LSC_1海深模型Fig. 12 BAT_LSC_1 bathymetry model

图选项

2.3 模型精度评价分析

引入ETOPO1和DTU18海深模型评估BAT_ LSC_1海深模型效能；同时采用未考虑地壳均衡补偿的重力异常导纳函数理论[36]恢复目标海区海深模型，建立的海深模型称为BAT_AF_1模型，研究海区各海深模型数值统计结果见表 4。

表 4 海深模型统计结果Tab. 4 Statistical results of bathymetry modelsm

表选项

表 4中各海深模型统计结果显示，研究海区海深最深约为5500 m，其中BAT_LSC_1最大海深处海深达到了6000 m；各模型海深平均值为4700 m左右；海深标准差均处于750 m左右，说明本文利用重力异常数据反演海深在试验海区与目前国际发布的海深模型表现相似。进一步评价本文海深模型质量，以事先准备的外部实测海深数据作为外部检核参考，采用双线性插值方法将海深模型值内插到检核点处，然后与检核点处海深作差并统计差值结果，各海深模型检核精度统计结果见表 5，表中相对精度为检核均方差与检核点平均海深绝对值之比。

表 5 海深模型检核差值统计结果Tab. 5 Statistical check results of bathymetry modelm

表选项

表 5检核差值统计结果显示，ETOPO1海深模型检核差值最大值达到了900 m，其余海深模型检核差值最值为500 m。就差值标准差和均方差统计指标而言，DTU18在研究海区检核精度最高，约为百米；BAT_LSC_1和BAT_AF_1检核精度优于ETOPO1模型，相较于ETOPO1，BAT_LSC_1检核精度提高了38%左右；BAT_AF_1检核精度提升约37%；BAT_LSC_1检核精度与BAT_AF_1相当，检核精度略高于BAT_AF_1。4种海深模型与检核点相关系数均达到了0.94以上，其中BAT_LSC_1和BAT_AF_1与检核点相关系数均高于0.97。ETOPO1检核相对精度为5.73%，BAT_LSC_1和BAT_AF_1相对精度均优于4%。综合海深模型检核差值统计分析结果，验证了本文方法构建海深模型技术方案可行，即利用统计方法构造地形-重力互协方差函数方法有效。进一步比对各海深模型检核效果，统计不同检核差值(差值取绝对值)范围内检核点所占比例，结果如图 13所示。

图 13 不同检核差值比例Fig. 13 The ratio of checkpoints in different ranges of difference

图选项

对比不同海深模型不同检核差值范围内检核点比例可以看出，DTU18海深模型明显优于其他3种海深模型。差值小于120 m左右范围时，ETOPO1模型检核点数量多于BAT_LSC_1和BAT_AF_1；检核差值大于120 m左右时，相同差值情况下ETOPO1模型检核点数量明显减少，检核点数量明显小于BAT_LSC_1和BAT_AF_1；检核差值范围大概120~280 m时，BAT_AF_1检核点数量多于BAT_LSC_1，之后BAT_LSC_1检核点数量略优于BAT_AF_1。经以上分析可知，试验海区利用最小二乘配置方法与利用导纳函数(未考虑地壳均衡补偿)方法构建的海深模型检核结果大体相当，说明了本文方法的实际可操作性。

2.4 陌生海区海底地形模型构建

据以上研讨分析结果可知，采用本文方法可有效构建目标海区海底地形数值模型，根据所建立的海底地形模型可较为清晰地反映目标海区地形地貌特征。最小二乘配置技术恢复海底地形模型关键核心是构造地形-重力协方差矩阵和重力异常自协方差矩阵。依统计分析理论，协方差矩阵反映了目标海区相关数据信息分布的内在规律，因此若某一陌生海区数据分布与已知海区数据分布具有相似特性，则可将反映已知海区数据分布内在特性的协方差函数迁移、推广、应用于陌生海区，进而实现以重力异常为输入的无船测数据支撑环境下，运用最小二乘配置方法反演海底地形目的。

选择西太平洋某5°×5°(157°E—162°E，19°N— 24°N)海域作为陌生备选海区，以上述1°×1°(155.2°E—156.2°E，22.5°E—23.5°E)试验海区为已知海区(以下称为区域1)。将陌生备选海区S&S V32.1重力异常作为已知数据(图 14)，相应的数据统计结果见表 6，通过比较、评判已知海区和陌生海区重力异常数据相似性程度，搜寻出备选海区中与已知海区数据特征最为相近区域。本文借鉴影像纹理特征匹配相关知识，首先将已知海区和备选海区重力异常作灰度化处理，得到像素值范围在[0, 255]的标准灰度影像，影像结果如图 15所示。

图 14 重力异常Fig. 14 Gravity anomaly

图选项

表 6 备选海区重力异常Tab. 6 Gravity anomaly in alternative sea areamGal

表选项

图 15 标准灰度影像Fig. 15 Standard grayscale image

图选项

为确定陌生备选海区中与已知海区(区域1)相似区域准确位置，同时兼顾影像匹配效率，本文分粗匹配和精匹配两个步骤环节在备选海区完成相似性区域搜索。首先，以已知海区重力异常灰度影像和陌生备选海区重力异常灰度影像为输入，根据SAD粗匹配算法模型，在备选海区中搜寻与已知海区最为相似的5幅子图，寻出的子图位置如图 16中小图的方框区域。然后，按照灰度共生矩阵(GLCM)和灰度-梯度共生矩阵(GGCM)生成原则，构造已知影像和5幅子图的GLCM和GGCM，考虑到计算消耗，将灰度值压缩为16个灰度级(具体构造过程见前文描述)，分别选择GLCM特征值中和平均、和方差、和熵和反差，GGCM特征值中灰度熵、灰度平均和灰度均方差，以及Hu矩(7个不变矩)一共14个特征值作为影像特征向量构成元素；以设计的特征向量为参考，分别计算5幅子图与已知影像间的欧式距离，与已知影像距离最小的子图即为最相似海区，精匹配结果如图 16(a)中蓝色方框所示区域(以下称为区域2)。另外，采用以上粗匹配和精匹配方法获得与已知海域相似度最低(与已知海域欧式距离最远)区域如图 16(b)中蓝色方框所示(称为区域3)。

图 16 影像匹配结果Fig. 16 Image matching results

图选项

从S&S V32.1重力异常模型中提取出区域2和区域3卫星测高重力异常数据，相应的统计结果见表 7，区域2和区域3重力异常平均值分别为4.13 mGal和25.36 mGal，标准差f分别为43.20 mGal和58.96 mGal。笔者收集到区域2实测水深数据分布如图 17中红色圆点所示，共计265个；S&S V19.1海深模型中通过实测水深数据约束生成的模型点共309个(认为S&S V19.1中实测水深值约束生成的模型点为实际水深测量值)，分布情况如图 17所示，相应的统计结果见表 7。可以看出，S&S V19.1中实测水深数据相较于笔者收集的船测数据更为丰富，因此将S&S V19.1实测水深数据作为评估海深模型精度的外部检核参考。区域3 S&S V19.1实测水深数据统计结果见表 7。

表 7 区域2、3数据Tab. 7 Data in area 2, 3

表选项

图 17 S&S V19.1实际船测水深数据和NGDC船测水深数据分布比较Fig. 17 Comparison of shipborne bathymetry data between S&S V19.1 and NGDC

图选项

据以上研究分析，因区域1和区域2海底地形地貌相似程度高，以区域2反演波段(15~160 km) S&S V32.1重力异常为输入，将区域1事先解算的地形-重力协方差函数、重力-地形协方差函数和重力异常协方差函数迁移于区域2，应用最小二乘配置技术可获得区域2反演波段海底地形；而后以S&S V19.1作为区域2先验海深模型，并施以滤波操作，从而可补充非反演波段(波长小于15 km和波长大于160 km)海底地形信息；叠加区域2所有波段地形信息，进而可实现陌生海域(区域2)海底地形模型构建。采用以上海底地形构建策略，建立的区域2海深模型(称为BAT_LSC_2)效果如图 18所示。

图 18 BAT_LSC_2模型Fig. 18 BAT_LSC_2 model

图选项

将区域1相关协方差函数应用于区域2，采用最小二乘配置方法恢复的区域2海底地形模型(BAT_LSC_2模型)可明显反映区域2的海底地形地貌特征，区域2海区海底分布着5座大型海底高山。区域2海底地形的呈现状态定性表明：以海洋表面重力异常为相似性参考，将已知海区解算的协方差函数迁移、推广并应用于相似未知海区，采用最小二乘配置方法恢复得到的海底地形可较清晰地表达未知海区海底地形特征。为进一步评估该方法的适用性，依据重力异常导纳函数反演理论，建立区域2海底地形模型(下文简称BAT_AF_2模型)。以事先设置的外部船测检核点为参考，各海深模型的检核统计结果见表 8。

表 8 海深模型检核结果Tab. 8 Check results of bathymetry model

表选项

表 8统计结果表明，BAT_LSC_2检核统计结果指标优于其余模型，如以检核均方差为例，BAT_LSC_2检核均方差为254.07 m，而BAT_AF_2、ETOPO1和DTU18分别为386.62、454.87和308.66 m。BAT_LSC_2与检核点相关系数最高，达到了0.970 6。4种海深模型中，ETOPO1模型检核相对精度超过10%，其余模型相对精度均在10%以内，其中BAT_LSC_2模型相对精度为6.16%。基于以上简要分析，进一步定量阐明了本文将已知海区相关数据统计计算结果迁移、推广并应用于相似陌生海区方法可行，构建的海底地形模型具有较强表现能力。比较各海深模型不同检核差值范围内检核点数量占检核点总数情况，结果如图 19所示。

图 19 不同检核差值范围比例Fig. 19 The ratio of checkpoints in different ranges of difference

图选项

图 19各海深模型不同检核差值范围检核点数量占比曲线比对可以发现，DTU18模型在差值小于400 m左右范围表现优于BAT_LSC_2，而BAT_LSC_2表现整体优于BAT_AF_2和ETOPO1；BAT_AF_2和ETOPO1模型表现相当。经以上分析可知，就模型检核统计各指标和不同检核差值范围检核点比例而言，利用最小二乘配置方法构建的区域2海底地形模型明显好于基于重力异常导纳函数构建的海底地形模型，进一步说明了本文解决无船测水深数据支撑环境下海底地形模型构建具有可行性。

类似地，将区域1解算的地形-重力协方差函数、重力-地形协方差函数和重力异常协方差函数迁移于区域3，应用最小二乘配置技术亦可反演获得陌生区域3海底地形信息，称为BAT_LSC_3模型。同理，采用导纳函数方法获得的区域3海底地形模型命名为BAT_AF_3模型。BAT_LSC_3模型和BAT_AF_3模型统计结果及其外部水深检核统计结果分别见表 9和表 10。

表 9 海深模型Tab. 9 Bathymetry modelsm

表选项

表 10 海深模型检核结果Tab. 10 Check results of bathymetry model

表选项

表 9显示BAT_LSC_3模型和BAT_AF_3模型整体表现基本一致。由表 10可知，BAT_LSC_3模型和BAT_AF_3模型检核均方差分别为347.03 m和359.56 m，BAT_LSC_3模型精度略优于BAT_AF_3模型。相较于BAT_LSC_2，BAT_LSC_3质量明显下降：①BAT_LSC_3系统误差明显增大，BAT_LSC_2检核差值平均值为46.08 m，BAT_LSC_3检核差值平均值达到了182.73 m；②BAT_LSC_2检核精度明显优于导纳函数结果，而BAT_LSC_3与导纳函数海底地形结果比较接近。不同检核差值范围内BAT_LSC_3模型和BAT_AF_3模型检核点数量占检核点总数情况如图 20所示。

图 20 不同检核差值范围比例Fig. 20 The ratio of checkpoints in different ranges of difference

图选项

对比图 20和图 19可知，BAT_LSC_3模型与导纳函数反演结果较为接近，检核差值小于250 m左右，BAT_AF_3模型表现明显优于BAT_ LSC_3；而BAT_AF_2模型则在整个检核差值范围表现均好于导纳函数结果。笔者认为原因是：使用最小二乘配置方法构建区域2和区域3海底地形模型时，均使用了区域1地形-重力协方差函数、重力-地形协方差函数和重力异常协方差函数进行解算，而区域2和区域3分别为与区域1相似度最高和最低区域，从而区域2协方差函数与区域1更为接近，进而区域2解算结果质量更高。另外值得说明的是，本文利用重力异常恢复海底地形的最小截止波段为15 km，即奈奎斯特(Nyquist)波长为15 km；则依据频率域采样定理可知，恢复的海底地形分辨率最高为30 km，而本文构建的海底地形模型在波长小于15 km部分利用了S&S V32.1模型补充，因此构建的海底地形模型最终分辨率依赖S&S V32.1模型，海底地形网格大小为1.8 km左右。

3 结论

本文尝试探索基于已知海区数据统计信息恢复未知陌生海区海底地形的手段途径。首先，利用已知海区可用的稀疏船测水深数据和卫星测高重力异常数据，统计建立了数据已知海区(区域1)地形-重力互协方差函数和重力自协方差函数，依据最小二乘配置方法恢复了区域1的海底地形；然后，以卫星测高海面重力异常作为相似性参考，采用影像粗匹配和精匹配算法模型，寻出了与已知海区重力异常最为相似的无船测水深数据海区(区域2)；最后，将区域1相关数据的协方差函数迁移、推广、应用于区域2，采用最小二乘配置方法构建了相应的海底地形模型，并评估了所构建的海底地形模型效能。经过定性和定量两个层面比对分析试验结果可以发现：本文方法反演的海底地形模型可有效描述研究海区海底地形地貌特点，反演的研究海区海深模型整体呈现效果好于根据重力异常导纳函数构建的海深模型，同时本文方法可巧妙回避导纳函数涉及的地球物理参数获取难题。试验区域1中，BAT_LSC_1模型检核精度相较于ETOPO1提升了约38%；试验区域2中，无论检核点差值统计结果，还是不同检核差值范围检核点比例情况，BAT_LSC_2表现均明显优于BAT_AF_2和ETOPO1。

作者简介

第一作者简介：范雕(1991-), 男, 博士, 讲师, 研究方向为物理大地测量和空间化海洋测绘工程。E-mail: fandiao2311@mails.jlu.edu.cn

通信作者：李新星, E-mail：minibad@126.com

初审：张艳玲复审：宋启凡
终审：金 君

资讯