, ]( B" ?% u; G# T
本文内容来源于《测绘学报》2024年第2期（审图号GS京(2024)0297号）
利用扰动重力数据反演海底地形的高斯曲面函数估计方法
' h7 a7 ?8 a" g 翟振和
, 孙中苗, 管斌, 马健, 李端
* `8 H) x r% C) t/ S; J. `, N& g 地理信息工程国家重点实验室, 北京 100029
9 _( n$ a4 [% {4 O( W
基金项目：国家自然科学基金（42174001；41774018）；地理信息工程国家重点实验室基金（SKLGIE2023-22-5）
; {# [) n. |+ h/ {6 I
摘要：利用2017年1月—2020年12月的SARAL卫星测高数据反演得到南海海域1×1格网扰动重力数据，通过与船载重力数据比较，精度达到5.5 mGal。提出了利用重力数据及高斯曲面函数求解区域特征参数进而估计海底地形模型的方法，利用SARAL卫星反演获得的重力数据在南海海域开展了计算试验。结果表明，利用ETOPO-1先验模型，在10×10格网一组划分条件下估计得到的1×1海底地形精度相对于先验模型提高了约10 m；利用DTU18先验模型，在9×9格网一组划分条件下估计得到的1×1海底地形精度相对于先验模型提高了约9 m。计算结果表明，利用卫星测高反演得到的重力数据，通过高斯曲面函数解算获得的5个特征参数，可以在一定程度上代表相应区域海底地形的曲面特征，进而可在不依赖外部实测数据条件下对先验海底模型进行精化。对于曲面估计而言，格网划分越小，曲面函数就越能反映区域变化特征。因此，对于未来卫星测高技术发展而言，更高分辨率的重力场探测技术有望继续提升海底地形细节的反演能力。
$ @; Z0 p6 P& D5 v& L! B) x) i
关键词：扰动重力 海底地形 SARAL卫星 高斯曲面估计 ETOPO-1 DTU18
翟振和, 孙中苗, 管斌, 等. 利用扰动重力数据反演海底地形的高斯曲面函数估计方法[J]. 测绘学报，2024，53(2)：231-238. DOI: 10.11947/j.AGCS.2024.20230204ZHAI Zhenhe, SUN Zhongmiao, GUAN Bin, et al. An estimation method of seabed topography based on Gauss surface function using ocean gravity data[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(2): 231-238. DOI: 10.11947/j.AGCS.2024.20230204 阅读全文：http://xb.chinasmp.com/article/2024/1001-1595/20240203.htm
9 G+ c0 ]' j/ S# Z# @9 r/ j 引 言* P- d# |9 }$ A8 N# k( B
海底地形信息对于船艇航行、地质资源勘查及地球科学研究等具有重要意义，目前，海底地形信息的获取主要采用船载多波束测量、机载激光测深、多光谱反演及卫星测高反演等技术[1-6]。船载多波束测量精度最高，已经成为海域高精度海底地形测绘的主要手段，但受限于测量平台和海洋客观环境条件，船载测量技术很难达到全球均匀覆盖。海洋卫星测高通过双星串飞、合成孔径高度计等手段可在2.5年左右时间获取全球高达1′分辨率的重力数据[7-9]，进而可反演获得一定分辨率且均匀分布的海底地形，现已成为海底地形三维模型构建的重要手段之一。我国自主研制的新一代海洋测高卫星于2023年3月发射，有望进一步提升宽阔海域海底地形信息的精细程度。目前，各国研究机构广泛应用包括船载测深、卫星测高等数据构建海底地形模型，美国国家环境信息中心(National Centers for Environmental Information，NCEI)、丹麦科技大学(Technical University of Denmark，DTU)及斯克里普斯海洋研究所(Scripps Institution of Oceanography，SIO)等多个研究机构相继发布了ETOPO-1、DTU18、SIO V19.1等全球海底地形格网模型[10-11]。国内学者基于卫星测高数据和船载水深数据也获得了局部或全球海域海底地形，如BAT_WHU2020[12]。在海底地形反演方法上，目前广泛利用卫星测高获得的重力异常数据通过重力地质法、解析法、频域法及最小二乘配置法等反演海底地形[13-25]，上述方法一般采用线性化后取1阶项进行计算，文献[26]顾及高阶项影响构建了非线性2次相关模型，利用该模型可在已知少量船载水深/重力数据条件下反演得到一定精度的海底地形。尽管利用卫星测高数据反演海底地形取得了很大进展，但考虑到重力场与海底地形之间的复杂函数关系，仍有必要加大研究并拓展反演思路。本文基于重力场与海底地形的随机相关特点，提出通过构建曲面特征的思路在先验模型基础上进行推估得到新的海底地形，为进一步拓展卫星测高数据的应用价值提供借鉴。
+ B$ i; b0 \8 a. e6 m6 l
/ [6 \& T& Z3 S6 d 1 考虑奇异影响的海洋扰动重力反演模型
2 S( ^: s1 U1 e4 W2 l' |! k
利用卫星测高数据反演重力场的一般方法是在垂线偏差基础上反演得到重力异常，考虑到扰动重力在大地边值问题上的潜在优势，本文采用扰动重力精确快速反演方法[27]计算得到格网点的扰动重力δgp
+ S! V. e: v) ^7 ?& v$ C' P
(1)
+ t2 d+ M: f6 ?" P! L: U
式中，γ表示正常重力；ξ、η分别表示垂线偏差的子午分量和卯酉分量；α表示计算点到观测点的方位角，核函数K(ψpq)的导数表示为
: I9 M8 i0 D5 ?1 X6 x# C" Y 
(2)
) ^6 f' P$ Z4 k( v+ L
当球面角距ψ趋近0时，式(2)会出现奇异问题，为了避免奇异，仿照文献[2]的思路在忽略二阶及以上项影响条件下推导得到扰动重力在奇异点处的计算公式为
7 k1 T' I" [: d7 _ 
(3)
) ?: T3 }# k; r6 j/ k }6 ~4 U q% `) ~ 式中，ξy、ηx分别表示垂线偏差子午分量在纬度方向的梯度、垂线偏差卯酉分量在经度方向的梯度；s0表示格网半径的大小。
: f6 Z: `: j# L) O, l
2 海域重力数据反演海底地形的高斯曲面函数估计方法
: {* \3 K* I) P% | 在已有学者研究中可以发现，在局部区域重力扰动/异常与地形起伏之间有较强相关性，为此本文提出这样一种假设，即在小区域重力数据的曲面特征与对应区域海底地形的曲面特征具有较强一致性，如果已知精确的重力场特征信息，就有可能恢复出相应的海底地形特征，并以此可估计出小区域的海底地形。为了对曲面变化特征进行定量分析，本文采用高斯曲面函数进行描述。二维高斯曲面函数可以用来描述两个维度的数据分布情况，在图像滤波、边缘检测等方面发挥着重要作用，可以生成高斯卷积核，用于图像处理中的高斯滤波，可有效去除高斯噪声，在机器学习核数据挖掘中也能发挥分类、聚类和回归等作用，从而更加准确地预测数据。高斯曲面函数具体表现为5个特征参数即区域峰值、区域峰值对应的经度、纬度，以及区域数据在经度和纬度方向上的变化幅度。
& h* }/ R! Z0 t4 j' q! E! H- g 假设卫星测高反演或其他手段获得的区域格网重力数据为δg(i)(本文设定格网大小为1′×1′)，i=1，2, …, n。将数据按照从北至南、从西至东划分为由m×m个1′×1′格网组成的小区域(m>3)。在每个小区域，利用已知重力数据按照高斯曲面函数模型(式(4))及最小二乘原理求解得到5个特征参数GT、GxT、GyT、CGx、CGy
1 t# k) F$ i. I# }% c: C3 [ t* b6 ^ e
(4)
$ Y6 D! _( T9 G! k4 X5 t( H1 a1 `
式中，GT表示局部区域内重力场(如重力异常、扰动重力)的峰值；GxT、GyT表示峰值对应的经度、纬度信息；CGx、CGy表示经度、纬度方向重力场变化的幅度；xi、yi表示每个重力场变量对应的经纬度信息。
' M$ r8 M0 w+ C, l
引入先验海底地形模型Hs(i)，i=1，2, …, n，格网大小为1′×1′。利用先验海底地形数据按照高斯曲面函数模型(式(5))及最小二乘原理求解得到5个特征参数HT、HxT、HyT、CHx、CHy
5 R; L v3 U" z* ^! y$ u! S, R 
(5)
" @7 g% ^1 ? `8 b+ k 式中，HT表示局部区域内海底地形的峰值；HxT、HyT表示峰值对应的经度、纬度信息；CHx、CHy表示经度、纬度方向海底地形变化的幅度；xi、yi表示每个海底地形变量对应的经纬度信息。
' E& P: W& U' z 从已有陆地重力、高程及海洋重力、水深数据(含格网数据和测线数据)分析来看，两种数据在变化的峰值位置是基本一致的，因此可将GxT、GyT表示待估计海底地形峰值对应的经度、纬度，将CGx、CGy赋以一定比例因子(式(6)、式(7))，而后作为待估计海底地形在小区域的变化幅度
4 ~# S# f4 R" N% v( ~$ L 
(6)
) V0 q% M- M. ^ 
(7)
9 a. u& ~7 _& K0 ]
将待估计海底地形的峰值HTnew在(HT-1000，HT+1000)范围内循环，将scale在(-1/10，1/10)范围内循环，上述循环为双重循环同时进行，循环次序对结果没有影响。而后按照式(8)进行求解
! g8 H4 @: Q3 o
(8)
7 ?+ K2 \ y: {, d- S 当小区域估计海底地形Hest－i和先验海底地形数据Hs(i)之间的均方差小于一定阈值时，循环结束，输出小区域估计海底地形。
4 U: ]" u8 f2 ?8 w7 |4 X 重复以上步骤，可将所有小区域计算完毕，最终输出得到完整区域的海底地形，整个估计方法的计算流程如图 1所示。
) G& E+ t& E0 D( R9 b
& J; k& z- K% E S; [2 k. j* A4 w 5 K4 H$ x0 P$ M& c4 W7 _
图 1 海底地形反演计算流程Fig. 1 Flowchart of seabed topography inversion
+ F1 p5 u) X. @4 j) b 图选项
& [2 k$ ]. C$ Y$ U4 y
! v7 L+ ]8 C8 {
% |8 U$ u2 }! Y: q9 b. i+ N
3 计算试验
, G6 e; N0 o4 }3 Y8 s! `8 u6 ]
3.1 海域扰动重力反演
# A1 k0 J" q" P% u: f( P& E
试验数据采用法国空间局发布的SARAL卫星GDR-F版数据，时间跨度为2017年1月—2020年12月，该时间段内SARAL卫星处于漂移轨道，因此数据较为密集。反演使用数据的区域范围为0°N—20°N，105°E—125°E，输出重力产品的区域范围为111.4°N—112.9°N，10°E—12.4°E，在该区域有高精度船载重力、水深测量数据(来自自然资源部)，分辨率达到1′×1′。为了分析SARAL卫星数据的海面分布密度，首先将官方发布的1 Hz海面高数据在重力场计算的局部海域进行画图显示，分别如图 2、图 3所示。
3 R2 t% S) z! J5 ^% m$ O( C5 v0 m 9 u2 R) _& Y. N4 d
N: }* z3 i6 O# E3 b2 `2 N 
图 2 SARAL卫星1 Hz海面高数据在1′×1′格网内的分布(蓝色圆点代表测量值)Fig. 2 Distribution in 1′×1′ grid of 1 Hz sea surface height data from the SARAL satellite (blue dots represent measurements). G E( j, M8 B3 ]" V) G/ w
图选项
. q5 d' }" ^+ i& p r. x
' K4 f9 T* U6 j3 r
# O# n6 C3 p6 J$ G, L" S; c 7 h; v& x) P0 Q/ A& Y
6 e' }; n. t6 Q1 l
图 3 SARAL卫星1 Hz海面高数据在3′×3′格网内的分布(蓝色圆点代表测量值)Fig. 3 Distribution in 3′×3′ grid of 1 Hz sea surface height data from the SARAL satellite (blue dots represent measurements)7 N3 Y6 |0 k) ~4 G: z+ `
图选项
4 v! Q7 G8 |' ?2 y w( h$ }
6 R8 L0 H- u- |- W" K
( ]5 h7 _8 z5 Y* H 由图 2、图 3可知，SARAL卫星1 Hz海面高数据在1′×1′格网内的分布很不均匀，不能完全达到1′×1′分辨率要求，总体上测量点格网分辨率优于3′×3′。实际重力场反演分别按照1′×1′、3′×3′格网大小进行计算，这样做的原因主要考虑以下3个原因：一是为了与已有1分格网船载重力数据进行比较，方便给出SARAL卫星反演重力场的精度；二是同时为了方便利用国际已经发布的1分格网的ETOPO-1、DTU18模型；三是考虑到本文方法需要在小区域利用已知格网观测量进行求解，而理论上格网划分越细，效果越好，同时考虑到目前试验区已知的船载水深格网模型为1′×1′格网，因此，采用1′×1′格网进行计算有利于对方法进行验证。计算中为避免远区项影响，采用移去恢复方法，先验模型采用2160阶EGM2008模型，按照式(1)—式(3)计算得到南海区域扰动重力如图 4所示。
" w+ }5 m# u- d! g
% }( s+ w4 G& F9 k/ D9 O6 e " ?* k7 }1 F5 | }3 _
图 4 SARAL卫星测高数据反演获得的1′格网扰动重力数据Fig. 4 The 1′ grid disturbing gravity data obtained from SARAL satellite altimetry data( k( C) h! ?; ~: |
图选项
% K9 B; @; H2 I s. G6 j' k
r) m" r$ \5 O: m% U# J5 x
2 o- }1 g+ i) Z, V 经与该区域船载重力数据比较，利用SARAL卫星测高数据及EGM2008模型得到的1′×1′格网扰动重力数据精度达到5.5 mGal，3′×3′格网扰动重力数据精度达到4.8 mGal(表 1)，符合现有卫星测高反演重力场的正常水平。
0 m0 n0 ~4 ]4 O. i 表 1 SARAL卫星测高数据反演扰动重力精度Tab. 1 Estimation accuracy of disturbing gravity inversion from SARAL satellite
5 F! P& ~/ I% M, f8 t
+ w. u7 {; h: y; P9 g" K: P* L
3.2 区域海底地形反演
% j( e' s* g8 |6 |5 V+ q2 M 试验区域和重力数据反演范围一致，具体为10°N—12.4°N，111.4°E—112.9°E，试验区平均海深3828 m，最大深度4316 m，最小深度1549 m，变化幅度562 m。采用的重力数据为上一节得到的1′×1′格网扰动重力数据。先验海底地形模型分别采用美国国家环境信息中心发布的1′×1′格网ETOPO-1水深模型(与船载水深比较结果见表 2)和丹麦科技大学发布的1′×1′格网DTU18海深模型(与船载水深比较结果见表 2)，船载水深检核格网数据在自然资源部2018—2019年的南海实测数据(测线间距约2 km)基础上格网化得到，分辨率达到1′。
7 M# M& N; i. ] 表 2 先验水深模型的精度统计Tab. 2 Estimation accuracy of seabed topography based on prior bathymetric model
4 P) {% C4 m( ]
4 J% O1 G8 b5 F9 R% I' i
由于曲面函数需要估计5个参数，这就使得重力观测格网数据有最小解算量，过少的观测量估计的参数准确度偏低，过多的观测量又会降低参数的适应范围，从而降低对高频信息的恢复。为了考察不同格网大小划分对计算结果的影响，将区域格网划分设定为3×3格网(即每个小区域有9个1′×1′格网)、4×4格网、5×5格网、6×6格网、9×9格网、10×10网格和12×12网格，不同格网划分条件下的计算结果见表 3，其中基于10×10网格划分反演得到的海底地形三维图形如图 5所示。
" t, G% |' F' C0 o g: i
表 3 基于ETOPO-1先验水深模型的海底地形估计精度Tab. 3 Estimation accuracy of seabed topography based on ETOPO-1 prior bathymetric model
0 S3 Y6 d4 m, w
* \5 J. c/ C: g2 a% |$ c/ O
, Z! P0 U4 c, ~4 N- s- h$ d# |/ j " r3 r5 U: r, Y. P/ u3 |4 t; G
图 5 基于ETOPO-1先验水深模型估计得到的海底地形Fig. 5 Seabed topography estimated based on ETOPO-1 prior bathymetric model) d' K/ Z. Q2 _. b) Q5 @( P
图选项
9 V8 c: n+ @. L( z' b
- h+ W: w3 Q# v+ S" q; _5 d. I' X
- W% H5 G7 F y# r( Q% r# p
基于上述方法及DTU18先验水深模型得到的估计海底地形，与船载数据比较得到结果见表 4，其中基于9×9网格划分反演得到的海底地形三维图形如图 6所示。
; C; W# |9 q) ?8 f* f1 m8 L) J 表 4 基于DTU18先验水深模型的海底地形估计精度Tab. 4 Estimation accuracy of seabed topography based on DTU18 prior bathymetric model
/ @! N2 T- B3 _; U4 }* x' Z
# N3 K3 D2 A. V
, Z/ d, ? c+ o9 O0 G# n2 J 3 ~, d e9 i% O% }, m. U: E2 \
图 6 基于DTU18先验水深模型估计得到的海底地形Fig. 6 Seabed topography estimated based on DTU18 prior bathymetric model
9 @! |& S& k, Y5 g h6 z$ L 图选项
- O8 ` y) x N& f5 F$ s% ?9 y
" e. G2 y" p& G' Y2 \ b 0 W! k6 X" F, s+ W5 p: q9 v+ I
从以上结果可以看出，虽然理论上大于5个观测量的重力数据都可以进行高斯曲面估计，但实际上，受制于重力数据本身分辨率及高斯曲面函数的非线性特征，其5个参数的求解受到数据量及数据分布影响还是比较大的，在网格数据划分上，数量过少或过多都不能取得好的效果。此外，如果海底地形出现复杂的变化趋势如非光滑、非连续界面，则高斯曲面函数受限于重力观测数据的分辨率难以重构出复杂的海底地形特征。综合分析看，对于ETOPO-1先验模型，10×10格网划分取得的效果最好，估计得到的海底地形精度相对于ETOPO-1提高了约10 m。对于DTU18模型，9×9格网划分取得的效果最好，估计得到的海底地形精度相对于DTU18提高了约9 m。综合来看，对于1′×1′格网海底地形反演而言，按照9×9格网划分进行分区求解，在不依赖外部海深数据支持下可以取得较优结果。
: y3 a/ P0 {/ F6 L) z5 m( E
4 结语
2 a9 ^% V c) H h/ p$ k
本文从重力数据和海深数据的随机相关特征出发，探索了利用卫星测高反演重力数据及曲面特征函数估计海底地形的方法，得到以下有益结论。
9 m- x( q4 [% a, n/ ?* W c5 d, X
(1) SARAL卫星测高1 Hz数据在南海区域的总体格网密度优于3′×3′，但不能满足1′×1′格网分辨率的要求，利用船载重力数据进行评估，其反演3′×3′格网重力场精度达到4.8 mGal，反演的1′×1′格网重力场精度优于5.5 mGal。
6 G. a1 ]4 A# `' y j, B' L1 K (2) ETOPO-1、DTU18都是国际上发布的比较权威的海深模型，其模型构建本身就采用了卫星反演水深和船载测量水深等多种来源的数据，以上述模型为先验模型，利用卫星测高获得的重力数据在9×9格网划分形式下通过高斯曲面函数获得特征参数进而估计海底地形，其精度较先验海深模型有约10 m量级的提升，这说明利用最新重力观测数据在不依赖外部实测水深数据条件下仍然可以对ETOPO-1、DTU18这些已有模型进行精化，也说明利用高斯曲面函数特征参数进行估计的思路是可行的。
, n/ V& u0 r7 }
(3) 总体而言，本文方法更适合在那些船载测量数据无法有效覆盖或覆盖不均匀的海域，在这些区域可利用高质量的重力数据并结合先验信息进行优化提升。实际计算中，高斯曲面函数的估计精度其实和观测数据的密集程度及精度有关，如果观测数据分辨率很高，则高斯函数的适应性也会很强，曲面函数就越能反映海底地形区域变化的细部特征(高频信息)。因此，重力场越精细，理论上所能反映海底地形的细部特征就越明显。对于未来以海洋测绘为任务目标的卫星测高技术发展而言，在提高精度的同时更应提高海面高的真实物理分辨率，这样才有可能利用高分辨率重力场去提升未测海域精细海底地形的反演能力。
作者简介
' o$ c2 K; N G- { 第一作者简介：翟振和(1980—), 男, 博士, 副研究员, 研究方向为空间大地测量与垂直基准。E-mail: zhaizhenhe1980@163.com
1 ?- ~, s, l+ H6 P4 V' a* f 初审：张 琳9 b1 l) W: s3 Z( A, Y& B
复审：宋启凡
& _6 O1 @. t( [5 V
终审：金 君
+ J7 F" H, C9 ]3 @7 V( Y
资讯
: s G1 x. \! j9 I; n) R
