2. 武汉大学 测绘学院,湖北 武汉,430079
2. School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China
1 引 言
全球海底地形模型是进行物理海洋学、地球物理学和海底地质学等地球科学研究的必要数据。但是,经过几十年的海洋测深数据积累,船测海深的覆盖面依然很稀疏,并且早期测深数据的精度也不高,根据船测海深构建的海底地形模型,其分辨率和精度均较低,例如GEBCO模型[1, 2, 3]。卫星测高技术的发展,为海底地形模型的构建提供了全新技术手段。文献[4]首先证明可以由测高数据推算海底地形。利用卫星测高技术,学者们计算获得了高精度海面重力异常数据[5, 6],可以将其应用于海底地形反演。接下来的研究表明[7, 8, 9, 10, 11, 12, 13, 14, 15],分辨率约为15~20km、相对精度约5%的海底地形模型,可由船测海深和卫星测高数据联合构建。当前,几乎所有相关研究机构发布的海底地形模型均是利用测高重力异常反演计算的,在特定波段内依赖于重力异常,一般不能直接将其应用于海底地壳均衡的三维导纳分析。文献[16]首次提出利用垂直重力梯度异常数据反演海底地形,并指出垂直梯度异常数据可以削弱海底地形反演时地壳均衡现象、沉积层等的影响,而且反演结果将不直接依赖于重力异常。但是,由于垂直梯度异常数据本身信噪比很低,文献[16]提出的空域反演法并不实用。文献[17]联合垂直梯度异常数据和ETOPO2海底地形模型,在频域内,进行了海底地形反演的试验性研究,但效果并不理想。文献[18]采用模拟数据,研究了利用垂直梯度异常反演海底地形时地壳均衡现象等因素的影响。根据垂直重力梯度异常反演海底地形的较优成果,目前在国内外尚未见发表。
本文采用斯克里普斯海洋研究所(Scripps Institute of Oceanography,SIO)发布的垂直重力梯度异常数据,联合来自美国国家地球物理数据中心(National Geophysics Data Center,NGDC)的船测海深数据,反演计算了全球75°S—70°N范围内的1′×1′海底地形模型。以印度洋南部和西北太平洋地区反演结果为例,通过与船测海深比较,分析了反演结果的精度,探讨了高次项和地壳均衡现象对反演结果的影响,并分析了联合重力异常、垂直重力梯度异常和船测海深反演海底地形的可行性。
船测海深来自NGDC,共3871航次,其分布情况如图 1所示。在大洋区域,尤其是南部大洋区域,船测海深的覆盖非常稀疏,且其中有很大一部分数据是卫星导航技术出现之前观测的,存在较大误差,已按文献[1]给出的处理方法对其进行了处理。1′×1′垂直重力梯度异常数据来自SIO,版本V20.1[19],是根据卫星测高海面倾斜数据计算获得的[20] 。卫星测高数据的海面轨迹间距一般为4~8km,在全球范围内,数据分布较船测数据密集得多,可以用于反演中短波段(<200km)的海底地形。
2.2 方 法首先,考察海底地形起伏与海面垂直重力梯度异常的关系。
根据地壳均衡理论,产生海面垂直重力梯度异常的最大异常物质源,是海底地形起伏及其地壳均衡补偿物质。根据岩石圈挠曲均衡模型,在频域内,由海底起伏引起的洋壳挠曲量为[21]
式中,k=2π/λ为波数;λ为波长;R(k)为地壳挠曲量的傅里叶变换;H(k)为海底地形起伏的傅里叶变换;ρm、ρc和ρw分别为地幔、地壳和海水密度;岩石圈挠曲响应函数为 式中,g为平均重力加速度;D为岩石圈抗挠刚度,D=ETe3/[12(1-υ2)],E为杨氏模量,一般取1011N/m2,Te为岩石圈有效弹性厚度,υ为泊松比,一般取0.25。根据文献[22]公式(4)和傅里叶变换的导数定理,同时顾及海底地形起伏及其均衡补偿物质的影响,直接给出频域内垂直重力梯度异常与海底地形之间的函数关系为
式中,ΔGz(k)为垂直梯度异常的傅里叶变换;G为万有引力常数;d为平均水深;t为平均洋壳厚度;F表示傅里叶变换;h(x)为空域内海底地形起伏;其他符号意义同前。式(3)即为由垂直梯度异常反演海底地形的理论基础。一般情况下,海底地壳密度、厚度、岩石圈有效弹性厚度等参数未知,且地壳均衡一般影响长波部分海底地形的构建精度,而长波部分海底地形可以由船测海深构建,因而,在较短波长(<200km)部分,式(3)中均衡影响部分可以忽略;再者,由于海深的向上延拓效应,式(3)中的高次项影响一般也可忽略,在海山地区,海底地形与海面垂直重力梯度异常的相干性如图 2所示。
从图 2看,在20~200km波长范围内,海底地形与垂直梯度异常数据的相干性很强,即表明垂直梯度异常信号在该波段内主要来源于海底地形。长波部分(>200km)受地壳均衡作用的影响,短波部分(<20km)受向上延拓效应、观测噪声的影响,想干性均减弱。因此,在20~200km波段内部分,可用垂直梯度异常反演海底地形,且两者之间的响应函数关系可简化为
当垂直重力梯度异常已知时,海底地形起伏的反演计算公式为 由式(5)可知,海底地形与垂直梯度异常经向下延拓等处理后的结果呈线性关系。但是,由于海底地质构造复杂,且海底地壳密度未知,其线性比例系数一般不由式(5)中的理论参数直接计算,而是根据计算区域的船测海深和垂直梯度异常数据进行计算[7],具体计算过程,将在下文中结合实例予以说明。本文海底地形反演结果将由两部分构成:一是由稀疏分布的船测海深构建的长波海底地形模型;二是由垂直重力梯度异常反演计算的20~200km(这一截断波长是参考重力异常反演文献,并结合笔者的试验计算确定的)波段海底地形。数据处理流程如图 3所示。
以印度洋南部区域(40°E—100°E,30°S—60°S)为例,海底地形的构建过程主要分为3步(本文数据处理过程中采用了FFT方法,因而实际数据处理都是以2°×2°分块进行的):
(1) 长波海底地形和反演波段内垂直梯度异常模型的构建。由船测海深(少量无船测数据的区块加入了ETOPO1数据)经格网化后[23],进行200km低通滤波处理,即可获得长波海底地形hlong(x),如图 4(a)所示。
垂直梯度异常数据经20~200km带通滤波、向下延拓等处理后,可获得反演波段内的垂直梯度异常模型Δgz(x),如图 4(b)所示。
(2) 反演波段内海底地形与垂直梯度异常之比的计算。根据式(5),在反演波段内,Δgz(x)与海底地形呈线性关系,但因海底各地区的地质构造状况不同,比例系数需根据船测海深计算。
计算过程可分3步进行:
第1步,根据长波海底地形(图 4(a)),采用双线性内插方法,内插计算船测点上的参考海深href(x′),并从观测值h(x′)中减去此参考值,得到船测点上的残差海深hres(x′),即
式中,x′表示船测点坐标。第2步,根据反演波段内垂直梯度异常模型(图 4(b)),内插船测点的垂直梯度异常Δgz(x′)。
第3步,计算船测点上的海底地形与垂直梯度异常之比
将S(x′)格网化,获得反演波段内的海底地形与垂直梯度之间格网化比值S(x),如图 4(c)所示。从图 4(c)看,在海盆区域,海底地形较平坦,海底地形与垂直梯度异常的关系较弱,因而比例系数接近零;而在海山、洋中脊等海底地形起伏较大的地区,比例系数则较大。(3) 计算海底地形模型。海底地形模型反演结果如图 4(d)所示,它由两部分构成,一是根据船测海深构建的长波海底地形hlong(x)(图 4(a));二是反演波段内的海底地形,该部分是比例系数(图 4(c))与垂直梯度异常模型(图 4(b))之积,即
式中,hpredict(x)表示海底地形反演结果。 3 计算结果及精度采用第2节所述的数据源和数据处理方法,计算了全球75°S—70°N范围内1′×1′海底地形模型,如图 5所示。计算过程中,三分之二的船测海深参与了计算,另三分之一的数据用于对计算结果进行精度检核,图 5中反演计算范围以外的海底、陆地及水深浅于100m区域的地形数据均取自V15.1模型。
海底地形模型构建的精度水平,与船测海深的稀疏程度有关,南部大洋船测海深数据较稀疏,因此本反演结果的精度在南部大洋相对北部大洋略低,在图 4(d)所示的印度洋南部地区,船测点上模型内插值与船测值之差,统计参数见表 1。
地形模型 | 最小值/m | 最大值/m | 平均值/m | 标准差/m | 均方根/m | 相对精度/(%) |
反演值 | -3200.405 | 2993.629 | 2.746 | 188.845 | 188.864 | 5.036 |
GEBCO | -3774.510 | 2396.339 | 18.969 | 239.004 | 239.755 | 6.373 |
DTU10 | -3827.382 | 2183.604 | 13.693 | 223.300 | 223.719 | 5.955 |
ETOPO1 | -2803.194 | 2602.072 | -7.073 | 214.190 | 214.306 | 5.712 |
V15.1 | -2754.816 | 2602.902 | -1.010 | 180.226 | 180.228 | 4.806 |
表 1中,参与比较的海底地形模型包括ETOPO1、V15.1、来自国际海道测量组织的GEBCO模型和来自丹麦科技大学的DTU10模型。从表中统计数据看,在印度洋南部,本文反演结果的精度与来自SIO的海底地形模型V15.1相当,优于GEBCO、DTU10和ETOPO1模型精度。检核点上残差频率分布直方图如图 6所示。从图中可知,残差主要分布在±100m范围内,残差值在±250m范围内的比例达到92%,仅3.9%的检核点残差绝对值大于400m。个别点的残差绝对值在1000m以上。这主要有两方面的原因:一是船测海深存在导航误差,二是在海底地形变化较剧烈的地区(如海槽边缘),采用双线性内插法计算的模型内插值可能存在偏差。
在北部大洋,船测海深数据较密集,海底地形反演计算的精度将更高,以西北太平洋为例,在0°—30°N,150°E—210°E范围内,本文反演结果如图 7所示,它与船测海深之差,统计见表 2,检核点水深残差频率分布直方图,如图 8所示。
地形模型 | 最小值/m | 最大值/m | 平均值/m | 标准差/m | 均方根/m | 相对精度/(%) |
反演值 | -2543.361 | 1991.111 | -0.146 | 132.908 | 132.908 | 2.740 |
GEBCO | -4893.587 | 4145.417 | -18.227 | 401.357 | 401.770 | 8.275 |
DTU10 | -3389.516 | 3506.325 | -17.423 | 295.131 | 295.645 | 6.085 |
ETOPO1 | -2566.856 | 2943.186 | 5.175 | 177.887 | 177.962 | 3.668 |
V15.1 | -3664.840 | 3277.342 | 1.223 | 158.513 | 158.517 | 3.268 |
从表 2看,在西北太平洋地区,本文海底地形反演的相对精度达到2.740%,优于ETOPO1、DTU10和GEBCO和V15.1模型。从图 8看,西北太平洋地区的残差分布较印度洋南部(图 6)更密集,残差绝对值在250m以内的点数占总检核点的95.8%,仅1.8%的残差绝对值大于400m,说明该地区的反演精度较印度洋南部更高。
4 讨 论 4.1 高次项及地壳均衡的影响根据反演得到的海底地形模型,采用CRUST2.0提供的洋壳厚度和密度参数,正演计算海底地形及其均衡补偿物质产生的垂直重力梯度异常(仅顾及(3)式中的n=1项),再在垂直梯度原始观测值中减去此计算值,将差值处理后,乘以海底地形构建过程中的比例系数,即可近似认为是高次项带来的最大影响。在印度洋南部,由高次项计算得到的海底地形均方根为92.404m,而反演结果的均方根为3867.681m,亦即高次项影响的相对量级最大仅为约2.4%。在西北太平洋地区,由高次项计算得到的海底地形均方根为95.981m,而海底地形反演结果的均方根为5017.344m,相对量级约为1.9%。可见,在海底地形构建过程中,可以略去式(3)中的高次项。
顾及地壳均衡影响时,洋壳密度结构参数同样来自CRUST2.0。在印度洋南部,取海底岩石圈有效弹性厚度Te=7km进行计算,计算结果与不顾及均衡影响时反演结果之差值的均方根为102.594m,相对量级为约2.7%。以船测海深为检核,标准差为201.873m,略低于不顾及均衡时的计算结果。在西北太平洋地区,取Te=5km进行海底地形反演计算,计算结果精度为199.168m,与不顾及均衡时反演结果之差值的均方根为102.684m,相对量级为约2.0%。可见,地壳均衡现象对海底地形反演精度的影响很小。在全球海底地形反演计算过程中,由于CRUST2.0的分辨率较低,且岩石圈有效弹性厚度未知,若顾及地壳均衡影响,一方面会增加计算量,另一方面不会提高计算精度,因此,采用本文所述的数据处理方法构建全球海底地形模型时,不必顾及地壳均衡现象的影响。
4.2 联合重力异常和垂直梯度反演海底地形垂直梯度数据具有抑制长波信号误差和放大短波信号的双重效果。在短波部分,本文通过数据计算试验,选定20km截断的高斯低通滤波器对数据进行滤波处理,在抑制高频噪声的同时,降低了反演模型的分辨率和精度。为考察垂直梯度异常数据对长波部分和短波部分反演结果的不同影响,在印度洋南部和西北太平洋地区,笔者联合重力异常和垂直梯度异常数据,进行了海底地形的反演计算。计算过程中,长波海底地形依然由船测海深格网化获得,100~200km波段范围内的海底地形由垂直梯度异常数据反演计算,小于100km波长的海底地形由重力异常数据反演计算。
在印度洋南部地区,联合重力异常和垂直梯度反演的海底地形精度为184.795m;在西北太平洋地区,反演精度有所提高,达到111.687m,较V15.1的精度(见表 2)提高了约29.5%。在两个区域内,采用重力异常反演、垂直梯度异常反演和联合两种数据的反演结果与船测海深之差的标准差见表 3。
地区 |
|
西北太平洋 | 149.221 | 132.908 | 111.687 |
印度洋南部 | 188.657 | 188.845 | 184.795 |
从表 3看,在西北太平洋地区,联合两种数据源反演计算的海底地形精度最高,根据垂直梯度异常反演计算精度优于根据重力异常反演计算精度。西北太平洋地区的海山分布密集(图 8),地形起伏较大,且覆盖了厚度约0.1~0.5km的沉积层[24]。表 3中的结果说明:在海底地形反演过程中,垂直梯度异常数据在100~200km波段内的表现较重力异常好,而波长小于100km部分则采用重力异常反演更优。文献[7]指出,海底沉积层质量异常是重力异常反演海底地形时的较大误差来源,而文献[16]认为垂直重力梯度异常数据可以抑制这一质量异常导致的海底地形反演时的不利影响。从表 3的结果看: 在100~200km波段内,采用垂直梯度异常数据反演时,相对较好地抑制了沉积层等中短波长质量异常的影响,提高了计算精度;在短波部分(<100km),由于垂直梯度异常数据信噪比较低,采用20km截断波长对数据进行低通滤波处理,在抑制了高频噪声的同时,也可能使海底地形的分辨率和精度略有降低。
在印度洋南部,3种反演方法获得的海底地形精度较为一致。这是因为,本地区海底地形起伏较小,海山顶部离水面较深(意味着无论对重力异常还是垂直梯度来说,短波地形产生的信号均较弱),沉积层主要分布在南部地形较平坦的海盆地区(地形平坦,易于通过稀疏的船测数据得到精度较高测长波海底地形),洋中脊地区无沉积层覆盖[24]。
4.3 短波端的截断及海底地形的分辨率问题当前,由SIO、DTU等机构发布的海底地形模型,一般是联合测高重力异常和船测海深数据构建的,其模型分辨率一般为1′×1′,但是,这并非海底地形数据的真实分辨率。海底地形的真实分辨率与测高数据密度、信噪比有关,同时也受海深延拓效应的影响。为抑制短波段噪声的影响,文献[7]采用维纳滤波器对重力异常数据进行处理
式中,A为滤波器设计参数;d为平均水深。文献[14]在根据其最新重力异常数据构建海底地形模型V15.1的过程中,取A=3891km4,此时,平均水深分别为2km、4km、6km时对应的分辨率(使W(k)=0.5)分别为12km、16km和20km,如图 9所示,这可以认为是V15.1模型的分辨率。本文在海底地形反演过程中,对垂直重力梯度异常数据进行了20km低通滤波处理,因而可以认为反演结果的分辨率为20km。从第3节的计算结果看,20km的低通滤波不会使海底地形反演精度急剧降低。
5 结 论本文联合船测海深和垂直重力梯度异常数据,构建了全球海底地形模型,给出了频域内海底地形起伏及其均衡补偿与海面垂直重力梯度异常之间的函数关系,详细阐述了数据处理流程,反演计算了全球1′×1′海底地形模型。以船测海深为检核参考,在印度洋南部区域,反演结果的精度与V15.1一致,相对精度约为5%;在西北太平洋地区,反演精度略优于V15.1模型。通过对计算结果精度的考察和讨论,本文可以得出如下结论:
(1) 联合垂直重力梯度异常数据和船测海深,可以获得精度较高的海底地形模型,本文计算的模型精度在印度洋南部地区与V15.1相当,在西北太平洋地区略优于V15.1。
(2) 海底地形反演过程中,地壳均衡现象和公式(3)中高次项的影响可以忽略。
(3) 垂直重力梯度异常数据在100~200km波段范围内表现较优,如结合重力异常和垂直梯度异常数据,海底地形的反演精度可进一步提高。
(4) 本文对垂直梯度异常数据采用的20km低通滤波处理,是通过数值计算试验得到的较合理值,在抑制高频噪声的同时,不会使反演精度急剧降低。
[1] | SMITH W H F. On the Accuracy of Digital Bathymetric Data[J]. Journal of Geophysical Research, 1993, 98(86): 9591-9603. |
[2] | SMITH W H F, SANDWELL D T. Conventional Bathymetry, Bathymetry from Space, and Geodetic Altimetry[J]. Oceanography, 2004, 17(1): 8-23. |
[3] | GOODWILLIE A M. User Guide to the GEBCO One Minute Grid[EB/OL].2008[2013-03-22].http://www.gebco.net. |
[4] | DIXON T H, NARAGHI M, MCNUTT M K, et al. Bathymetric Predictions from SEASAT Altimeter Data[J]. Journal of Geophysical Research, 1983, 88(C3):1563-1571. |
[5] | HUANG Motao, ZHAI Guojun, GUAN Zheng, et al. On the Recovery of Gravity Anomalies from Altimeter Data[J]. Acta Geodaetica et Cartographica Sinica, 2001, 30(2): 179-184. (黄谟涛,翟国君, 管铮, 等. 利用卫星测高数据反演海洋重力异常研究[J]. 测绘学报, 2001, 30(2): 179-184.) |
[6] | LI Jiancheng, NING Jinsheng, CHEN Junyong, et al. Determination of Gravity Anomalies over the South China Sea by Combination of TOPEX/Poseidon, ERS2 and Geosat Altimeter Data[J]. Acta Geodaetica et Cartographica Sinica, 2001, 30(3): 197-202. (李建成,宁津生,陈俊勇,等. 联合TOPEX/Poseidon,ERS2和Geosat卫星测高资料确定中国近海重力异常[J]. 测绘学报,2001, 30(3): 197-202.) |
[7] | SMITH W H F, SANDWELL D T. Bathymetric Prediction from Dense Satellite Altimetry and Sparse Shipboard Bathymetry[J]. Journal of Geophysical Research: Solid Earth, 1994, 99(B11): 21803-21824. |
[8] | SMITH W H F, SANDWELL D T. Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings[J]. Science, 1997, 277(5334): 1956-1962. |
[9] | HWANG C W. A Bathymetric Model for the South China Sea from Satellite Altimetry and Depth Data[J]. Marine Geodesy, 1999, 22(1): 37-51. |
[10] | CALMANT S. MURIEL B N, CAZENAVE A. Global Seafloor Topography from a Least-square Inversion of Altimetry-based High-resolution Mean Sea Surface and Shipboard Soundings[J]. Geophysical Journal International, 2002, 151(2): 795-808. |
[11] | WANG Yong, XU Houze, ZHAN Jingang. High Resolution Bathymetry for China-sea and Adjacent Waters[J]. Chinese Science Bulletin, 2001, 46(11): 956-960. (王勇,许厚泽,詹金刚. 中国海及其邻近海域高分辨率海底地形[J]. 科学通报, 2001, 46(11): 956-960.) |
[12] | HUANG Motao, ZHAI Guojun, OUYANG Yongzhong, et al. The Recovery of Bathymetry from Altimeter Data[J]. Geomatics and Information Science of Wuhan University, 2002, 27(2): 133-137. (黄谟涛,翟国君,欧阳永忠,等. 利用卫星测高资料反演海底地形研究[J]. 武汉大学学报: 信息科学版, 2002, 27(2): 133-137.) |
[13] | LUO Jia, LI Jianchen, JANG Weiping. Bathymetry Prediction of South China Sea from Satellite Data[J]. Geomatics and Information Science of Wuhan University, 2002, 27(3): 256-260. (罗佳,李建成,姜卫平. 利用卫星资料研究中国南海海底地形[J]. 武汉大学学报: 信息科学版, 2002, 27(3): 256-260.) |
[14] | FANG Jian, ZHANG Chijun. 2′×2′ Sea Floor Bathymetry Prediction of China Sea and Its Vicinity[J]. Geomatics and Information Science of Wuhan University, 2003, 28(sup3): 38-40. (方剑,张赤军. 中国海及邻近海域2′×2′海底地形[J]. 武汉大学学报: 信息科学版,2003,28(增3): 38-40.) |
[15] | BECKER J J, SANDWELL D T, et al. Global Bathymetry and Elevation Data at 30 Arc Seconds Resolution: SRTM30_PLUS[J]. Marine Geodesy, 2009, 32(4): 355-371. |
[16] | WANG Y M. Predicting Bathymetry from the Earth’s Gravity Gradient Anomalies[J]. Marine Geodesy, 2000, 23(4): 251-258. |
[17] | WU Yunsun, CHAO Dingbo, LI Jiancheng, et al. Recovery of Ocean Depth Model of South China Sea from Altimetric Gravity Gradient Anomalies[J]. Geomatics and Information Science of Wuhan University, 2009, 34(12): 1423-1425. (吴云孙, 晁定波, 李建成, 等. 利用测高重力梯度异常反演中国南海海底地形[J]. 武汉大学学报:信息科学版, 2009, 34(12): 1423-1425.) |
[18] | HU Minzhang, LI Jiancheng, LI Dawei. Bathymetry Prediction from Vertical Gravity Gradient Anomalies[J]. Journal of Geodesy and Geodynamics, 2013, 32(5): 95-98. (胡敏章,李建成,李大炜. 利用垂直重力梯度异常反演海底地形[J]. 大地测量与地球动力学,2012,32(5): 95-98.) |
[19] | SANDWELL D T, SMITH W H F. Global Marine Gravity from Retracted Geosat and ERS-1 Altimetry: Ridge Segmentation Versus Spreading Rate[J]. Journal of Geophysical Research, 2009, 114(B01): 411-419. |
[20] | SANDWELL D T, SMITH W H F. Marine Gravity Anomaly from Geosat and ERS-1 Satellite Altimetry[J]. Journal Geophysical Research, 1997, 102(B5): 10039-10054. |
[21] | WATTS A B. Isostasy and Flexure of the Lithosphere[M]. Cambridge: Cambridge University Press, 2001. |
[22] | PARKER R L. The Rapid Calculation of Potential Anomalies[J]. Geophysical Journal of the Royal Astronomical Society, 1973, 31(4):447-455. |
[23] | WESSEL P. BECKER J M. Interpolation Using a Generalized Green’s Function for a Spherical Surface Spline in Tension[J]. Geophysical Journal International, 2008, 174(1): 21-28. |
[24] | DIVINS D L. Total Sediment Thickness of the World’s Oceans & Marginal Seas[EB/OL]. NOAA National Geophysical Data Center. Washington: NASA, 2003[2013-03-22]. http://www. ngdc.noaa.gov/mgg/sedthick/sedthick.html |