: L, `8 T+ ?" _/ o7 C+ x 【编者按】根据重力异常与海底地形的相关性，利用日本西之岛附近海域内的重力异常和船测水深，分别采用Sandwell&Smith方法（S&S，1994）和重力地质法（GGM）反演了该海域内1′×1′海底地形模型，并将反演的模型与格网化船测水深模型和ETOPO1模型进行了对比分析。在GGM中，使用了迭代法向下延拓确定密度差常数为5.7g/cm3。与S&S模型相比，GGM模型与格网化船测水深模型的差值的标准差达到了37.10m。本文发表在《海洋测绘》2016年第5期上，现编发给朋友们阅读了解。李倩倩，女，1989年出生，山东菏泽人，中国科学院测量与地球物理研究所，硕士研究生，主要从事测高重力场反演研究。 
9 B V8 y! ?9 e. r 一、引言 1 F& D- p' X7 V/ I( m3 C( Z4 c g5 R
海底地形作为海洋生物学、海洋化学、海洋动力学等海洋学科的基础，其信息的完备和精度具有重要的意义。在大地测量学的传统方法中，对海底地形的认识主要是通过船载回声测深实现的，但是由于传统的方法即耗时耗力，测量结果又分布有限、不均匀，所以由其测得的海底地形偏差很大。卫星测高技术的发展，为构建高精度的海底地形模型提供了一种新的、可靠的技术和方法。底地形在一定程度上基于重力异常与海底地形的相关性及由重力异常推估海深的研究近上世纪70年代陆续开展起来[1]。近几十年国内外已有学者对由重力值推算海深做过讨论研究[2-14]。Calmant对该领域的理论与方法进行了详细的总结[5]。其中经典方法以Smith&Sandwell方法（S&S）为代表，基于Parker公式[13]和Watts的挠曲均衡模型[14]，经过滤波、向下延拓等处理，在一定波长范围内确定重力异常与海深的线性关系，从而通过重力异常计算海深。另一种方法为重力地质法（GGM:Gravity-Geologic Method），海底洋壳和海水之间较小的密度差异变化使得重力地质法十分适宜于利用卫星测高重力异常反演海底地形的研究。现如今，GGM方法已经被成功地应用于南中国海、东日本海、格陵兰岛南部以及阿拉斯加南部等多处海域的地形反演[6-11]，在对格陵兰岛南部海域的海底地形反演中，将船测水深和反演结果进行比较，两者差异的标准差甚至达到±35.8m。 o# o% p1 r! n) i8 t" c% G
本文联合船测数据和测高重力异常数据，分别采用经典方法（S&S）和重力地质法（GGM）两种方法计算了西之岛附近海域（138°E～143°E，25°N～30°N）的1′×1′海底地形模型。在重力地质法中，本文采用向下延拓法确定了密度差常数为5.7g/cm3。此外，本文将上述两种方法反演的海底地形模型与船测水深和ETOPO1模型进行了对比分析。与S&S模型相比，GGM模型与格网化船测水深模型的差值的标准差达到了37.10m。
; m, x+ p5 t' c* f K 二、研究区域和数据
: _" ^* N, ~$ E# N, y1 i* ` 本文研究范围位于日本西之岛附近海域138°E～143°E，25°N～30°N之间，西之岛位于伊豆-小笠原-马里亚纳俯冲带，该俯冲带是太平洋板块向菲律宾海板块俯冲的边界。研究区域内区域地形复杂，含有大量海底火山、海沟、岛礁，海底起伏大，进行海底地形反演意义重大。本文采用的数据有以下2种：①船测海深数据，来自NGDC(National Geophysical Data Center)，本文采用的船测点共31897个，所有船测数据同时包含重力测量和水深测量。选取其中27341个作为控制点，剩余4556个为检核点，控制点与检核点不重复，见图1，其中黑色为控制点，白色为检核点。②卫星测高重力异常数据，来自SIO，UCSD(Scripps Institution of Oceanography, University of California,San Diego)，版本为V23.1，是2014年最新发布的，其联合了最新卫星测高资料，特别是加入了Jason-1末期大地测量任务周期数据和部分Cryosat-2资料，精度达到1～2mGal（1mGal=10-5m/s2)[4]，见图2。 
: T# ~: l) ? ^* O Y! g* w3 ? 图1 船测航迹分布  % A+ m* a3 i4 }6 O
图2 测高重力异常V23.1 # R4 h7 q5 u6 i) X: _. o+ k; _
三、GGM反演海底地形
. ?. l* j6 ` \; }. u2 D+ _) z ⒈ GGM基本原理 . t+ m8 A3 R" o. O& S0 ]7 D
GGM方法应用于海底地形反演的基本原理是：将重力异常划分为长波参考场和短波残差场两部分，其中长波部分由已知海深和重力异常的控制点计算，再在待求点上用重力异常减去长波重力异常得到短波分量，最后根据短波分量计算海深。我们设定，j1～jn为海域内的船测水深控制点，其重力异常和海深均为已知值；i1～in为待求点，其重力异常已知，海深待求。
' Y" [" @. K- `2 E# b* z9 Y: m 使用控制点jn的水深计算重力异常短波残余分量，公式如下：
- s9 b- }, p R: N, I+ b( h Δgshortjn＝2πGΔρ(Ejn－D) ⑴ ; f) M# V) Y4 a
式中，G为引力常数；Δρ为海水和岩石圈的密度差；Ejn为jn控制点的海深；D为参考海深，一般取所有船测控制点的最深值；Ejn和D的单位均为m。 ! C4 z4 m* P# D1 E! ?
将船测水深控制点观测重力异常减去其由公式⑴计算得到的短波部分得到控制点长波重力异常。通过适当的插值方法可求得待求点i1～in上的长波重力异常，再利用如下两式就得到i1～in点的短波重力异常和水深： : X( O# w( ]' ~* D$ N
Δgshorti＝Δgobsi－Δglongi ⑵ # \- P/ z% E& j, p
Ei＝D＋Δgshorti／2πGΔρ ⑶
2 d3 q% w( q4 \, _ 式中各符号意义同前。GGM方法计算海深时只涉及Δρ和D两个参数。D作为一个参考量，其取值对反演结果的影响很小。因此，海水和岩石圈的密度差Δρ的确定是该方法的关键。
/ s# Z" ^9 G$ _" d6 m. X ⒉ 向下延拓法确定密度差常数 : i0 M. E K* T" y: [
位场的向下延拓是一个病态问题，存在不稳定性，其解是不唯一的。本文采用了迭代法向下延拓[15]。相对于其他方法，迭代法原理简单，计算速度快。原理如下：设定平面ΓA与ΓB之间是无源空间，用u(x,y,o)表示ΓA平面（z＝0）的位，用U(kx,ky,0)＝F[u(x,y,o)]表示u(x,y,o)的傅里叶变换。ΓB平面（z＝h）上的位为： 
" U6 T. M+ ^ q) w- `9 k8 o 其中F-1是傅里叶变换，kx，ky是x，y方向的波数。式⑷是FFT的延拓公式。当h＞0时，上式是向上延拓；当h＜0时，上式是向下延拓，延拓不稳定。 " M8 O# }% U* @9 ^* X" d+ W0 ^/ M6 m
本文利用迭代法将海平面测高卫星重力异常向下延拓到西之岛附近海域最深约9789m处，迭代计算了50次。延拓面的重力异常与海平面的重力异常的比值等于延拓面的密度差与海平面的密度差（1.03g/cm3）之比[7]。根据这个原理，本文计算得到延拓面的密度差为5.7g/cm3，并选其进行GGM模型反演。
( N" j- F( X# w 四、S&S反演海底地形 6 T: g3 I1 S$ D0 s o
根据Parker公式，海底地形与重力异常之间的关系可以描述如下:
! v: d" ^* j" d% z H(K)＝Z-1(K)H(K)＝[1／(2πGΔρ] ×exp(2πkd)G(k) ⑸
|$ |/ K, f: C2 [$ S. \ 其中，G为地球引力常数；Δρ为岩石圈与海水的密度差；d是平均深度；k是径向频＝率（√kx2＋ky2)；K为(kx,ky)＝(1/λx,1/λy)，其中(kx,ky)和(λx,λy)分别表示x和y方向的频率和波长。G(K)和H(K)分别表示重力异常（Δg）和海底地形（h）的傅里叶变换；Z(k)表示响应函数。 * s: L% a& y& G8 q
由于上式中含有exp(2πkd)项，海洋深度反演实际上就演变成为一个向下解析延拓问题，这种问题的解是非唯一的，而且是不稳定的，因此，必须采取一些相应的规则化处理措施。其具体步骤参考文献[7]，本文不再赘述。 
' V. J: P3 t/ l H0 o4 r8 t (a)重力地质法反演 
3 p3 b3 j \" n1 n! Q, E$ K) I (b)经典方法反演 / Q3 Q* J9 w' n) a; p
图3 海底地形模型示意图 # m# J. F" E/ y8 E
表1 各模型的统计分析（单位：米） , w3 l) g2 h o7 l, ]6 ~5 c
Model 4 b% U- \+ `9 V, t D
Min ) {" T' G0 f9 \/ d( Z. ?
Max
& | g7 I! O5 [- | Mean 1 j( Z) e: I6 d' y e& i& o& ^
STD
5 A3 V' K* R9 d$ ` RMS ( Q3 j, K- ~& r8 x+ }4 k" {6 K A
Correlation coefficient(CC) $ H, D0 f7 G+ o N
GGM
9 z- P; {/ r9 j& a S&S
$ E% L! t6 [) @2 T% L1 g" a NGDC 6 W+ u$ q/ L" C5 s* Y2 S
ETOPO1 ' f# n5 v' r* O' f, Q& Z D
GGM ) k. {, V' ?, _
-9737.90
# @! ?" S' j% N2 s1 h 96.60 4 P9 c z2 r. T R* J
-3447.09 + D5 ?$ _6 b; H$ C
1355.94 q9 x2 @9 y. o: T' S
3704.18 2 H3 a' k4 Q C5 u* U0 ]
1
' C1 o# J5 `4 y) E7 X' |# ?& l$ p 0.9934 1 h) n0 }0 r3 v. J4 m
0.9996
- g' w0 ]; Q% F. I1 A+ x 0.9930
# {% f7 h* }) u, ~! P8 R: F S&S " v! d; N5 }% E! Y! ^9 v/ W
-9735.26
[& f( w& {+ c. V 547.75 : k+ N1 C4 p1 k0 o
-3463.87 6 D; M2 ?+ y$ g. s$ `
1353.22
# y) l$ i, X4 h# j 3718.82 ( V6 d5 ] r2 O7 e* B7 j. t
0.9936
9 E9 l C+ o _2 Z& z5 V; q. } 1
$ l8 S$ p" E( }; k2 o 0.9933
! j% s" ?' c# F% q% g 0.9873 # J3 ~1 c4 O+ n/ n$ e/ }
NGDC
& ^# g0 W0 s# D$ {- a -9769.25 ! [8 k2 W5 Q& ]4 I
-63.50 : ]9 p# X* q3 ?+ z& ~/ X/ ?8 l( O
-3458.98
, q5 T9 m- r9 p2 J$ V- H 1354.35
8 I7 U7 K) x* p 3714.68
" A. a& a4 l" @: s 0.9996 2 [" a: d y6 |8 J& z @
0.9933 3 F: }# O) ?& D, @7 T) o
1
( u5 W, m5 r' k' v4 Y 0.9926
. @' i& T. o9 C% m$ j) i( M+ ?' T ETOPO1
6 ^7 n3 @5 \6 \% f -9789.00
' h; Z' h; @" G' R 331.00
3 { C8 @& e4 {5 s. i -3478.10
' [( C) Z( ^+ M2 e 1375.44 : Y4 g/ ~- H% F+ V! ]
3740.19 8 H. f6 Q5 A, k/ ]
0.9930
, j w1 _ U- `/ V0 F5 {6 n6 n 0.9873
$ \1 a) I5 s9 ~4 i 0.9926 3 N" E1 c6 S/ G) s) `/ Y8 F' S
1 1 l3 g: b0 H. q+ w2 i! ]
图4给出了各模型与船测水深之间的差值图，其统计结果如表2所示。图4（a）所示的GGM模型与NGDC模型的差值图变化最小，对应表2中的最小STD（37.10m）。但图4（b）、（c）中的浅水区域（小笠原脊、Shichito-Iojima脊附近）差值很大，有的甚至达到±2000m，这是因为在浅水区域地形复杂，水深梯度变化剧烈。从表2中可以看出，相对于ETOPO1模型与NGDC模型差值的标准差168.50m，S&S模型的标准差为156.56m，这说明S&S模型与ETOPO1模型精度相当。综上所述，相对于S&S，GGM更适用于反演海底地形。  + \) s/ S" ~, N( X2 Q* v0 a6 O( B
(a)GGM-NGDC  & |8 [3 z! ~8 W" G) P- s
(b)S&S-NGDC  7 h, f' R7 M5 i. x# M% Q
(c)ETOPO1-NGDC
1 p. _0 Y$ O9 b; e! K" G" R W; G7 P 图4 各模型与船测水深之间的差值图 & g) K# m, v+ l
表2 各模型与船测水深之间的差值统计 单位：米 ! t& R: |+ T9 t0 l1 c% A, r9 B3 P% J
odel 3 B# O7 v9 ~; r/ L% s8 i
Min
9 e0 z& d! Y3 L* t* B. k Max : B" X9 Y# L5 t
Mean
( v. H3 R2 e" c5 B$ i- `6 Z3 R, w STD $ M3 c/ x1 y; a( M! \+ p% I
RMS
{* [4 z( E0 W. ` GM-NGDC 2 u' h8 ?, ]7 V# w- t/ K. n
-1164.05 2 z K0 T4 G7 @. i2 ~
893.80
" a2 `$ t( Q/ A, w/ A* Q 11.90 9 h" J) {" L' I/ W7 e
37.10 L# d$ l3 b, j$ Q! Q3 s8 M( l( |
38.97
6 S: H: @# L3 l& t GM-NGDC 2 G y2 z9 q+ D2 `% K% g7 l
-1164.05
! }* B* T8 O/ x- [$ `, _ j 893.80 2 \; F2 K4 z' [: e l) [3 |
11.90
9 d) K$ }( C- y. l5 ], t/ e% a 37.10 . `6 L; x3 f6 w
38.97 / S5 A2 e2 B5 E& F: H) }
ETOPO1-NGDC
! ?, t X% O- P9 }8 \( \) _ s5 k9 p -2221.16
' t! K- T3 O7 g$ x* ]' b. q9 Z 1842.03
8 m2 a8 T9 d4 L4 M -19.12 3 I% h2 H- @ X
168.50
/ h' |6 E* U1 u* A4 |* X 167.41
1 ^4 n+ \8 {% g 六、结论
; m2 Z1 h2 `7 _1 Z6 ? 本文根据重力异常与海底地形的相关性，联合船测数据和测高重力异常数据，分别采用经典方法（S&S）和重力地质法（GGM）两种方法计算了西之岛附近海域（138°E1～143°E，25°N～30°N）的1′×1′海底地形模型，并对其进行了精度评价，通过分析，得到以下结论： # m! _' O8 q* w! q
⑴采用重力地质法计算海底地形模型结果可靠，方法简便，能够有效填补稀疏海域水深测量的不足，不涉及复杂的地球物理参数求取；且与S&S模型相比，GGM模型在精度上具有明显的优势。
& O `9 M+ V- M. f0 X! \: K; x ⑵经典方法主要基于Parker公式，以船测数据为参考，直接计算海底地形与重力异常之间的线性比例系数，将重力异常转换成海底地形。该方法原理简单，计算过程中采用了FFT方法，速度较快。相比GGM模型，S&S模型精度较低，原因可能是在S&S方法的移去-恢复过程中，由于低于15km波长的短波分量中含有高频噪声，没有采取有效的措施将之恢复。
$ T$ O" K, C H' y) q ⑶GGM方法中，密度差Δρ的选择则是海底地形反演的关键。但是，由于受到各种地质因素的影响，密度差Δρ失去其本来的物理意义，作为一个经验参数参与计算。 
- F+ E; _, w/ _, v' y' Q 参考文献：
$ e( V8 A+ |9 d& {) g9 ^ [1]黄谟涛,翟国君,欧阳永忠,等.卫星测高技术应用研究回顾与展望[J].海洋测绘,2004,24(4):0065-0070. . S! j& G* }5 y) g" a! H- x1 x
[2]罗 佳,李建成,姜卫平.利用卫星资料研究中国南海海底地形[J].武汉大学学报:信息科学版,2002，27(3):256-259. , f' Z% I+ N2 l h
[3]Smith W H F，Sandwell D T.Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry[J].Journal of Geophysical Research,1994，99(B11):21803-21824.
1 x' I% Z8 l& w5 O0 S: |4 z [4]Sandwell D T,Muller D R，Smith W H F，et al.New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure[J].Science,2014,346(6205):65-67.
" ~' ^5 @! g6 t [5]Calmant S.Modeling bathymetry by inverting satellite altimetry data: a review[J].Marine Geophysical Research,1996,18:123-134. # c' p. O0 d* e r ]5 m* T
[6]欧阳明达,孙中苗,翟振和.基于重力地质法的中国南海海底地形反演[J].地球物理学报,2014,57(9): 2756-2765. & d" f$ J9 S/ q9 y$ w
[7]Kim K B,Hsiao Y S,Kim J W,et al.Bathymetry enhancement by altimetry-derived gravity anomalies in the East Sea(Sea of Japan)[J].Marine Geophysical Researches,2010,31(4):285-298.
$ ?5 p2 L5 C9 i6 y3 ]) y [8]胡敏章,李建成,金涛勇.应用重力地质法反演皇帝海山的海底地形[J].武汉大学学报（信息科学版）,2012,37(5):610-612,629.
! {6 s! j8 l, {& c- e/ _ [9]Jeong W K,Ralph R B,Bang Y L,et al. Altimetry-derived Gravity Predictions of Bathymetry by the Gravity Geologic method[J].Pure and Applied Geophysics.2011,168(5):815-826.
* o" j3 E( \$ b/ r B [10]Hsiao Y S,Kim J W,Kim K B,et al.Bathymetry estimation using the gravity-geologic method:An investigation of density contrast predicted by downward continuation method[J].Terr.Atmos.Ocean.Sci,2010,22(3):347-358. + E/ j) i& m( v3 U/ |3 m
[11]胡敏章,李建成,金涛勇.顾及局部地形改正的GGM海底地形反演[J].武汉大学学报（信息科学版）,2013,38(1):60-63.
) ~0 H5 l( g. J- v: Q% f4 ]% T [12]Hwang C.A Bathymetric Model for the South China Sea from Satellite Altimetry and Depth Data[J]. Marine Geodesy,1999,22:37-51.
5 {$ J+ z( y8 N, r$ I, { [13]Parker R L.The Rapid Calculation of Potential Anomalies[J]. Geophys.J.R.Astr.Soc.1972,31:447-455.
7 ?" @8 [2 O! f7 P8 V [14]Watts A B.An Analysis of Isostasy in the World’s Oceans 1:Hawaiian-Emperor Seamount Chain[J].Journal of Geophysical Research,1978,83:5989-6004.
N0 }! M l, Q, e8 D" ^* `8 y [15]徐世浙.迭代法与FFT法位场向下延拓效果的比较[J].地球物理学报,2007,50(1):285-289.
0 s5 j7 U* F" s4 A, R0 W/ c ■论文专区的文章均为在《海洋测绘》上刊发的论文，若其他公众平台转载，请备注论文作者，并说明文章来源，版权归《海洋测绘》所有。
9 t% W* F; Z+ A 相关阅读推荐：
9 a4 j; N- ]$ F; V" m0 C# x. E 论文专区▏海洋重力测量动态环境效应分析与补偿
( |7 v, \$ M* c$ k( n/ U; w 测绘百科▏重力异常与海洋重力基点 % M ~& v' V; h! v
论文专区▏多型航空重力仪同机测试及其数据分析 . k0 P" Q! R) H3 `
科技前沿▏动态重力仪技术的最新发展 + ^( j" ?( Y5 ?4 q2 `1 M; o- ?
海洋论坛▏南海珊瑚礁地貌的高分辨率遥感影像分类 0 c9 z, g. l0 v! e8 F" }3 Z" B, i
海洋论坛▏海洋遥感卫星发展历程与趋势展望 ( P- ^+ {$ n9 P9 X: ~; ]3 @+ l
论文专区┃河口海域地形变化检测技术探讨 0 Q" ]$ v& h$ V. d+ J/ u. @* V
测绘讲坛▏关于海底地形测量的认识 % Q7 W9 i) ^) ^- i8 |0 `( D
投稿邮箱：452218808@qq.com，请您在留言中标注为投稿，并提供个人简介及联系方式，谢谢！  o9 x3 ]5 x1 s7 w( W
1 B% q! S* `: w4 N
5 Q @3 o2 c4 e- U) c/ ?* @" ] y. K0 r! A9 _ }# {
2 X2 s6 k" f% x' S | 
