胡敏章 李建成 邢樂(lè)林 翟振和 歐陽(yáng)明達(dá)

1)中國(guó)地震局地震研究所(地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室)，武漢 430071

2)中國(guó)地震局地殼應(yīng)力研究所武漢科技創(chuàng)新基地，武漢 430071

3)武漢大學(xué)測(cè)繪學(xué)院，武漢 430079

4)武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079

5)西安測(cè)繪研究所，西安 710054

海底地形反演方法比較*

胡敏章1，2，3，4)李建成3)邢樂(lè)林1，2)翟振和5)歐陽(yáng)明達(dá)5)

1)中國(guó)地震局地震研究所(地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室)，武漢 430071

2)中國(guó)地震局地殼應(yīng)力研究所武漢科技創(chuàng)新基地，武漢 430071

3)武漢大學(xué)測(cè)繪學(xué)院，武漢 430079

4)武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079

5)西安測(cè)繪研究所，西安 710054

研究比較了海底地形反演的空域法、頻域法和最小二乘配置法，選擇3個(gè)2°×2°實(shí)驗(yàn)區(qū)域進(jìn)行海底地形實(shí)算。結(jié)果表明，空域法計(jì)算結(jié)果精度雖然較頻域法和最小二乘配置法略高，但計(jì)算速度非常慢，不適用于大區(qū)域海底地形模型的計(jì)算;頻域法由于可以采用FFT算法，計(jì)算速度很快，精度也較高;最小二乘配置法需要已知海底地形的詳細(xì)統(tǒng)計(jì)信息，這與海底地形反演目的相矛盾。本文還采用頻域法，聯(lián)合垂直重力梯度異常和船測(cè)海深數(shù)據(jù)反演了實(shí)驗(yàn)區(qū)域海底地形，驗(yàn)證了垂直重力梯度異常數(shù)據(jù)用于海底地形計(jì)算的可行性。

海底地形;重力異常;空域法;頻域法;最小二乘配置法

海底地形反演是衛(wèi)星測(cè)高的一大應(yīng)用領(lǐng)域，自Dixon等［1］首次采用Seasat沿軌大地水準(zhǔn)面數(shù)據(jù)計(jì)算海底地形以來(lái)，眾多學(xué)者都研究了由測(cè)高重力異常/大地水準(zhǔn)面反演海底地形的理論和方法，反演計(jì)算了大量局部或全球海底地形模型［2－15］。Calmant［11］對(duì)衛(wèi)星測(cè)高數(shù)據(jù)反演海底地形的理論和方法進(jìn)行了詳細(xì)總結(jié)。Wang［16］提出根據(jù)重力垂直梯度異常數(shù)據(jù)反演海底地形的方法，并認(rèn)為該方法計(jì)算得到的海底地形模型與重力異常相互獨(dú)立，能夠用于海底均衡狀況的三維導(dǎo)納分析。吳云孫等［8］對(duì)此進(jìn)行試驗(yàn)，但由于梯度數(shù)據(jù)短波部分噪聲很大，效果并不理想。胡敏章等［10］采用模擬數(shù)據(jù)，研究了根據(jù)垂直重力梯度異常反演海底地形的方法。

本文總結(jié)歸納現(xiàn)有海底地形反演方法，并通過(guò)實(shí)際數(shù)值計(jì)算，對(duì)比分析各方法的優(yōu)劣。分別根據(jù)重力異常和垂直重力梯度異常，采用頻域法反演了實(shí)驗(yàn)區(qū)海底地形，比較了兩種數(shù)據(jù)反演結(jié)果，說(shuō)明垂直梯度異常數(shù)據(jù)也可用于海底地形反演。

1 海底地形反演方法

1.1 空域法

如圖1所示，Q表示相對(duì)于參考海深的一個(gè)高為h的海底地形，a為地球平均半徑，Rref=a+d，d為參考海深，ψ為P與流動(dòng)點(diǎn)Q之間的球心角。海面上點(diǎn)P受Q處海底地形異常物質(zhì)產(chǎn)生的重力異常影響的計(jì)算公式為:

圖1 海底起伏的重力效應(yīng)Fig.1 Gravity effect of the seafloor undulation

式中，ΔgP表示P點(diǎn)周圍n個(gè)格網(wǎng)點(diǎn)的地形起伏在P 點(diǎn)產(chǎn)生的重力異常，ΔΩi、hi、ψi分別表示第 i個(gè)格網(wǎng)點(diǎn)的面元積分、地形起伏及其與P點(diǎn)的球心角距，面元積分近似為:

式中，θ為積分流動(dòng)點(diǎn)余緯，Δλ、Δφ分別為格網(wǎng)的經(jīng)度方向和緯度方向大小，φQ為格網(wǎng)中心緯度。

式中，B=R － acosψ，C=3 －4cosψ，D=1 －5cosψ +4cos2ψ，argsh是反雙曲正弦函數(shù)。海底地形與重力異常之間函數(shù)關(guān)系建立后，根據(jù)擬牛頓法，海底地形迭代計(jì)算模型為:

式中，hk為第k次迭代計(jì)算結(jié)果，h0為初始模型值，Chh為初始模型先驗(yàn)協(xié)方差矩陣，A為(1)式的線性化設(shè)計(jì)矩陣，E為觀測(cè)數(shù)據(jù)誤差矩陣，一般為對(duì)角矩陣，d為觀測(cè)數(shù)據(jù)向量，可以包括重力異常、大地水準(zhǔn)面、船測(cè)海深等多種數(shù)據(jù)，f(hk－1)表示對(duì)初始模型進(jìn)行非線性運(yùn)算，函數(shù)f將海深模型參數(shù)與觀測(cè)值聯(lián)系起來(lái)。

1.2 頻域法

根據(jù)Parker公式，在僅顧及一次項(xiàng)的情況下，海底地形起伏與重力異常之間的關(guān)系為:

式中，ΔG(k)為反演波段內(nèi)重力異常的傅里葉變換，ρc、ρw分別為洋殼和海水密度，dw為平均水深，H(k)為反演波段內(nèi)海底地形的傅里葉變換，其他符號(hào)意義同前。

采用頻域法聯(lián)合船測(cè)海深和重力異常數(shù)據(jù)反演計(jì)算海底地形模型的最終結(jié)果為:式中，h(x)為海底地形計(jì)算結(jié)果，hlong(x)為船測(cè)數(shù)據(jù)控制的長(zhǎng)波海底地形，S為反演波段內(nèi)海底地形與重力異常之比，Δgd(x)為反演波段內(nèi)向下延拓至平均水深處的重力異常。

1.3 最小二乘配置法

由重力異常計(jì)算海底地形的統(tǒng)計(jì)算法模型為［6］:

圖2 實(shí)驗(yàn)區(qū)域船測(cè)數(shù)據(jù)分布Fig.2 Distribution of the data of Ship Soundings in the Experimental Areas

此即為最小二乘配置公式［17］。Hwang［6］闡述了根據(jù)海底地形和重力異常的功率譜來(lái)估計(jì)Chg的方法，以確定配置法中重力異常至海底地形的轉(zhuǎn)換函數(shù)。式(9)變?yōu)?

形式上，式(10)與頻域法式(6)根據(jù)均衡原理得到的理論導(dǎo)納關(guān)系類似，但實(shí)際意義不同。Zcol(k)基于海底地形和重力異常的統(tǒng)計(jì)信息，將地形及重力異常均視為隨機(jī)信號(hào)，無(wú)地球物理意義。

2 算例分析

數(shù)值計(jì)算實(shí)驗(yàn)在3個(gè)2°×2°區(qū)域中進(jìn)行，區(qū)域一范圍28°～30°N，178°～180°W，共有 9 092 個(gè)船測(cè)海深觀測(cè)值(圖2(a));區(qū)域二范圍20°～22°N，156°～158°E，共有5 199個(gè)船測(cè)海深觀測(cè)值(圖2(b));區(qū)域三范圍38°～40°N，165°～167°E，共有 2 548個(gè)船測(cè)海深觀測(cè)值(圖2(c))。計(jì)算中長(zhǎng)波海底地形均以船測(cè)數(shù)據(jù)為基礎(chǔ)構(gòu)建，測(cè)高重力異常(或重力垂直梯度異常)僅用于計(jì)算波長(zhǎng)小于200 km部分的海底地形?？紤]到實(shí)驗(yàn)區(qū)域船測(cè)數(shù)據(jù)較少，將其中95%用于計(jì)算，另外的5%用于對(duì)計(jì)算精度的檢核。

2.1 空域法計(jì)算結(jié)果

以船測(cè)數(shù)據(jù)格網(wǎng)化模型作為初始模型，取海底地形先驗(yàn)?zāi)Ｐ驼`差和船測(cè)海深誤差為300 m，測(cè)高重力異常誤差為5×10－8ms－2，海底地形經(jīng)驗(yàn)協(xié)方差函數(shù)參數(shù)為0.25°時(shí)，根據(jù)式(5)僅進(jìn)行一次迭代計(jì)算，即可解得實(shí)驗(yàn)區(qū)域海底地形(圖3)。顧及計(jì)算速度和邊界效應(yīng)，每次計(jì)算時(shí)輸入數(shù)據(jù)范圍為0.5°×0.5°，取區(qū)域中心4點(diǎn)計(jì)算結(jié)果為有效值。2°×2°范圍內(nèi)共有14 400個(gè)數(shù)據(jù)點(diǎn)，共需進(jìn)行3 600次計(jì)算，以內(nèi)存為2 G的筆記本電腦計(jì)算，耗時(shí)約120 h，計(jì)算速度非常慢。

2.2 頻域法計(jì)算結(jié)果

根據(jù)頻域法，采用重力異常反演計(jì)算得到的實(shí)驗(yàn)區(qū)域海底地形如圖4所示。采用FFT算法，完成2°×2°范圍數(shù)據(jù)計(jì)算，僅耗時(shí)約1 min，計(jì)算效率大大優(yōu)于空域法。

同樣，利用垂直梯度異常數(shù)據(jù)替代重力異常數(shù)據(jù)，采用30 km短波段截?cái)嘁砸种聘哳l噪聲影響，即對(duì)數(shù)據(jù)進(jìn)行30～200 km帶通濾波后，反演計(jì)算的海底地形如圖5所示。

2.3 最小二乘配置法計(jì)算結(jié)果

以本文采用頻域法計(jì)算的初步結(jié)果(圖4)為起始模型，采用文獻(xiàn)［6］給出的功率譜密度函數(shù)法，計(jì)算反演波段內(nèi)的Zcol(k)。為方便應(yīng)用，以10階多項(xiàng)式對(duì)其進(jìn)行擬合(圖6)。

圖6中，紅色實(shí)線為根據(jù)數(shù)據(jù)計(jì)算得到的轉(zhuǎn)換函數(shù)，藍(lán)色虛線為10階多項(xiàng)式擬合結(jié)果。總體上，轉(zhuǎn)換函數(shù)呈帶通濾波形式，即通過(guò)對(duì)重力異常帶通濾波處理，可以將重力異常轉(zhuǎn)換成海底地形。根據(jù)圖6所示的轉(zhuǎn)換函數(shù)，在實(shí)驗(yàn)區(qū)域內(nèi)，計(jì)算結(jié)果如圖7所示。

2.4 精度檢核

圖3 空域法反演結(jié)果Fig.3 The bathymetry predicted by Space Domain Method

圖4 頻域法反演結(jié)果(基于重力異常)Fig.4 The bathymetry predicted by Frequency Domain Method(Based on gravity anomaly)

圖5 頻域法反演結(jié)果(基于垂直梯度異常)Fig.5 Bathymetry predicted by the frequency domain method(Based on vertical gravity gradient anomaly)

圖6 轉(zhuǎn)換函數(shù)Fig.6 Transform function

圖7 配置法計(jì)算結(jié)果Fig.7 Bathymetry predicted by least square collocation method

以船測(cè)數(shù)據(jù)為檢核標(biāo)準(zhǔn)，V15.1、ETOPO1、DTU10和GEBCO模型在區(qū)域一內(nèi)精度分別為115.8、234.8、342.1 和 435.7 m;在區(qū)域二內(nèi)精度分別為 65.3、302.7、360.8 和 442.8 m;在區(qū)域三內(nèi)精度分別為 27.8、137.0、188.3 和 180.5 m，4 個(gè)模型的精度依次呈下降趨勢(shì)，V15.1精度最高。3種方法計(jì)算結(jié)果與船測(cè)檢核數(shù)據(jù)之差的統(tǒng)計(jì)參數(shù)見(jiàn)表1。區(qū)域一、二、三的平均水深分別為約4 470、4 520和5 200 m，表1中，相對(duì)誤差為標(biāo)準(zhǔn)差與平均水深之比。V15.1來(lái)自斯克里普斯海洋研究所(SIO)，是采用頻域法，聯(lián)合船測(cè)海深和測(cè)高重力異常構(gòu)建的海底地形模型［2－3］;ETOPO1 是美國(guó)國(guó)家地球物理數(shù)據(jù)中心(NGDC)于2008年8月公布的全球地形模型［18］;GEBCO模型是在政府間海洋學(xué)委員會(huì)(IOC)和國(guó)際航道組織(IHO)主持下，主要由全球500 m等高距的數(shù)字化等深線數(shù)據(jù)以及船測(cè)海深數(shù)據(jù)，采用格網(wǎng)化方法獲得的［19］;DTU10模型是丹麥科技大學(xué)(DTU)2010年發(fā)布的全球海底地形模型，它基于測(cè)高重力異常和GEBCO模型構(gòu)建，20～120 km波段內(nèi)海底地形由重力異常反演計(jì)算，＞120 km波段海底地形及陸地地形均來(lái)自GEBCO模型，相關(guān)數(shù)據(jù)、文檔、軟件從 DTU 網(wǎng)站(www.space.dtu.dk)獲得。

從表1看，在區(qū)域一和區(qū)域二，空域法計(jì)算結(jié)果精度優(yōu)于ETOPO1、DTU10和GEBCO模型，低于V15.1模型，其他方法精度優(yōu)于DTU10和GEBCO模型，與ETOPO1精度一致;在區(qū)域三，空域法計(jì)算精度優(yōu)于DTU10和GEBCO模型，與ETOPO1精度一致，低于V15.1模型，其他方法計(jì)算結(jié)果精度均略低。

頻域法中根據(jù)重力垂直梯度異常反演計(jì)算的海底地形精度與根據(jù)重力異常反演計(jì)算的精度一致，說(shuō)明在實(shí)驗(yàn)區(qū)域內(nèi)采用30 km的短波段截?cái)啵谝种浦亓Υ怪碧荻犬惓?shù)據(jù)高頻噪聲的同時(shí)，沒(méi)有使計(jì)算精度顯著降低。

配置法計(jì)算精度與頻域法根據(jù)重力異常計(jì)算的精度較為一致，這是因?yàn)榕渲梅ㄞD(zhuǎn)換函數(shù)計(jì)算時(shí)采用的參考模型是頻域法結(jié)果。計(jì)算結(jié)果說(shuō)明，配置法計(jì)算精度受制于初始模型，為了獲得精度較高的計(jì)算結(jié)果，首先應(yīng)獲得計(jì)算區(qū)域內(nèi)較為可靠的海底地形統(tǒng)計(jì)信息。

3 結(jié)論

1)空域法雖然計(jì)算精度較高，但計(jì)算速度非常慢，不適用于大區(qū)域海底地形模型計(jì)算。

2)配置法的關(guān)鍵是計(jì)算轉(zhuǎn)換函數(shù)，這需要計(jì)算海底地形與重力異常之間的相關(guān)函數(shù)，因而需要預(yù)先已知海底地形較準(zhǔn)確的統(tǒng)計(jì)信息，這與海底地形反演相矛盾，實(shí)際中較少采用。

3)頻域法依據(jù)Parker公式，以船測(cè)數(shù)據(jù)為參考，直接計(jì)算海底地形與重力異常之間的線性比例系數(shù)，將重力異常轉(zhuǎn)換成海底地形。由于可以采用FFT方法，頻域法計(jì)算速度很快，精度也較好，是當(dāng)前常用的海底地形反演方法。V15.1和DTU10模型均采用了此方法。

4)垂直重力梯度異常數(shù)據(jù)也可用于聯(lián)合船測(cè)海深數(shù)據(jù)計(jì)算海底地形模型，且計(jì)算精度與重力異常反演一致。

1 Dixon T H.Bathymetric predictions from SEASAT altimeter data［J］.J Geophys Res，1983，88:1 563 － 1 571.

2 Smith W H F，Sandwell D T.Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry［J］.J Geophys Res，1994，99:21 803 －21 824.

3 Smith W H F，Sandwell D T.Global sea floor topography from satellite and ship depth soundings［J］.Science，1997，277:1 956－1 962.

4 Tscherning C C，Knudsen P，F(xiàn)orsberg R.First experiment with improvement of depth information using gravity anomalies in the mediterranean sea［R］.Thessaloniki:Department of Geodesy and Surveying，University of Thessaloniki，1994.

5 Calmant S，Berge-Nguyen M，Cazenave A.Global seafloor topography from a least-square inversion of altimetry-based high-resolution mean sea surface and shipboard soundings［J］.Geophys J Int，2002，151:795 －808.

6 Hwang C.A bathymetric model for the South China Sea from satellite altimetry and depth data ［J］.Marine Geodesy，1999，22:37 －51.

7 黃謨濤.利用衛(wèi)星測(cè)高資料反演海底地形研究［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2002，27(2):133－137.(Huang M T.The recovery of bathymetry from altimeter data［J］.Geomatics and Information Science of Wuhan University，2002，27(2):133 －137)

8 吳云孫.利用測(cè)高重力梯度異常反演中國(guó)南海海底地形［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2009，34(12):1 423 －1 425.(Wu Yunsun.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):1 423 －1 425)

9 李大煒.利用衛(wèi)星測(cè)高數(shù)據(jù)反演中國(guó)海及鄰近海域海底地形［J］.大地測(cè)量與地球動(dòng)力學(xué)，2009(4):70－73.(Li Dawei.Research on inversion of seafloor topography of China Sea and its adjacent areas from satellite altimetric data［J］.Journal of Geodesy and Geodynamics，2009(4):70 －73)

10 胡敏章.利用垂直重力梯度異常反演海底地形［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012(5):95－98.(Hu Minzhang.Bathymetry prediction from vertical gravity gradient anomalies［J］.Journal of Geodesy and Geodynamics，2012(5):95－98)

11 Calmant S.Modeling bathymetry by inverting satellite altimetry data:a review［J］.Marine Geophysical Research，1996，18:123 －134.

12 Watts A B.Isostasy and flexure of the lithosphere［M］.Cambridge University Press，2001.

13 Calmant S.Seamount topography by least-squares inversion of altimetric geoid heights and shipborne profiles of bathymetry and/or gravity anomalies［J］.Geophys J Int，1994，119:428－452.

14 Ramillien G，Wright C.Prediction seafloor topography of the New Zealand region:a nonlinear least square inversion of satellite altimetry data ［J］.J Geophys Res，2000，105(B7):16 577－16 590.

15 Parker R L.The rapid calculation of potential anomalies［J］.Geophys J R Astr Soc，1973，31:447 － 455.

16 Wang Y M.Predicting bathymetry from the earth’s gravity gradient anomalies［J］.Marine Geodesy，2000，23(4):251－258.

17 Moritz H.Advanced physical geodesy［M］.Wichmann Karlsruhe，1980.

18 Amante C，Eakins B W.ETOPO1 1 arc-minute global relief model:procedures，data sources and analysis［R］.NOAA Technical Memorandum NESDIS NGDC －24，2009.

19 Goodwillie A M.User guide to the GEBCO one minute grid［R］.2008.

COMPARATIVE ANALYSIS OF METHODS FOR BATHYMETRY PREDICTION

Hu Minzhang1，2，3，4)，Li Jiancheng3)，Xing Lelin1，2)，Zhai Zhenhe5)and Ouyang Mingda5)
1)Key Laboratory of Earthquake Geodesy，Institute of Seismology，CEA，Wuhan 430071
2)Wuhan Base of Institute of Crustal Dynamics，CEA，Wuhan 430071
3)School of Geodesy and Geomatics，Wuhan University，Wuhan 430079
4)Key Laboratory of Geospace Environment and Geodesy，Ministry of Education，Wuhan 430079
5)Xi’an Institute of Surveying and Mapping，Xi’an710054

Three methods of bathymetry prediction:Space Domain Method，F(xiàn)requency Domain Method and Least Square Collocation Method were compared by actual bathymetry calculating for three experimental areas.The results show that the accuracy of Space Domain Method is slightly better than the others，the calculation efficiency is much lower and not suitable for bathymetry predicting over a large area，however.Frequency Domain Method can be used to calculate bathymetry model quickly by FFT algorithm，and accuracy is acceptable.Detailed statistics information of the bathymetry is needed for bathymetry inversion with Least Square Collocation Method.This is contradicted with the purpose of bathymetry predicting.The bathymetry models for the experimental areas were calculated with Frequency Domain Method combined with the vertical gravity gradient anomaly and the data from ship sound-ings.The result indicates that it is possible to predict bathymetry wiyh vertical gravity gradient anomalies.

bathymetry;gravity anomaly;method in the space domain;method in the frequency domain;leastsquares collocation

P229.1

A

1671-5942(2014)05-0011-06

2013-11-09

中國(guó)地震局地震研究所所長(zhǎng)基金重點(diǎn)項(xiàng)目(IS201326125);國(guó)家測(cè)繪地理信息局測(cè)繪基礎(chǔ)研究基金項(xiàng)目(13－01－01);國(guó)家自然科學(xué)基金項(xiàng)目(41204019，41304003)。

胡敏章，男，1985年生，博士，助理研究員，研究方向?yàn)楹５椎匦闻c重力均衡。E－mail:huminzhang@126.com。