海底地形反演方法比較*

胡敏章 李建成 邢樂林 翟振和 歐陽明達

1)中國地震局地震研究所(地震大地測量重點實驗室)，武漢 430071

2)中國地震局地殼應力研究所武漢科技創新基地，武漢 430071

3)武漢大學測繪學院，武漢 430079

4)武漢大學地球空間環境與大地測量教育部重點實驗室，武漢 430079

5)西安測繪研究所，西安 710054

海底地形反演方法比較*

胡敏章1，2，3，4)李建成3)邢樂林1，2)翟振和5)歐陽明達5)

1)中國地震局地震研究所(地震大地測量重點實驗室)，武漢 430071

2)中國地震局地殼應力研究所武漢科技創新基地，武漢 430071

3)武漢大學測繪學院，武漢 430079

4)武漢大學地球空間環境與大地測量教育部重點實驗室，武漢 430079

5)西安測繪研究所，西安 710054

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

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

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

本文總結歸納現有海底地形反演方法，并通過實際數值計算，對比分析各方法的優劣。分別根據重力異常和垂直重力梯度異常，采用頻域法反演了實驗區海底地形，比較了兩種數據反演結果，說明垂直梯度異常數據也可用于海底地形反演。

1 海底地形反演方法

1.1 空域法

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

圖1 海底起伏的重力效應Fig.1 Gravity effect of the seafloor undulation

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

式中，θ為積分流動點余緯，Δλ、Δφ分別為格網的經度方向和緯度方向大小，φQ為格網中心緯度。

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

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

1.2 頻域法

根據Parker公式，在僅顧及一次項的情況下，海底地形起伏與重力異常之間的關系為:

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

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

1.3 最小二乘配置法

由重力異常計算海底地形的統計算法模型為［6］:

圖2 實驗區域船測數據分布Fig.2 Distribution of the data of Ship Soundings in the Experimental Areas

此即為最小二乘配置公式［17］。Hwang［6］闡述了根據海底地形和重力異常的功率譜來估計Chg的方法，以確定配置法中重力異常至海底地形的轉換函數。式(9)變為:

形式上，式(10)與頻域法式(6)根據均衡原理得到的理論導納關系類似，但實際意義不同。Zcol(k)基于海底地形和重力異常的統計信息，將地形及重力異常均視為隨機信號，無地球物理意義。

2 算例分析

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

2.1 空域法計算結果

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

2.2 頻域法計算結果

根據頻域法，采用重力異常反演計算得到的實驗區域海底地形如圖4所示。采用FFT算法，完成2°×2°范圍數據計算，僅耗時約1 min，計算效率大大優于空域法。

同樣，利用垂直梯度異常數據替代重力異常數據，采用30 km短波段截斷以抑制高頻噪聲影響，即對數據進行30～200 km帶通濾波后，反演計算的海底地形如圖5所示。

2.3 最小二乘配置法計算結果

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

圖6中，紅色實線為根據數據計算得到的轉換函數，藍色虛線為10階多項式擬合結果。總體上，轉換函數呈帶通濾波形式，即通過對重力異常帶通濾波處理，可以將重力異常轉換成海底地形。根據圖6所示的轉換函數，在實驗區域內，計算結果如圖7所示。

2.4 精度檢核

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

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

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

圖6 轉換函數Fig.6 Transform function

圖7 配置法計算結果Fig.7 Bathymetry predicted by least square collocation method

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

從表1看，在區域一和區域二，空域法計算結果精度優于ETOPO1、DTU10和GEBCO模型，低于V15.1模型，其他方法精度優于DTU10和GEBCO模型，與ETOPO1精度一致;在區域三，空域法計算精度優于DTU10和GEBCO模型，與ETOPO1精度一致，低于V15.1模型，其他方法計算結果精度均略低。

頻域法中根據重力垂直梯度異常反演計算的海底地形精度與根據重力異常反演計算的精度一致，說明在實驗區域內采用30 km的短波段截斷，在抑制重力垂直梯度異常數據高頻噪聲的同時，沒有使計算精度顯著降低。

配置法計算精度與頻域法根據重力異常計算的精度較為一致，這是因為配置法轉換函數計算時采用的參考模型是頻域法結果。計算結果說明，配置法計算精度受制于初始模型，為了獲得精度較高的計算結果，首先應獲得計算區域內較為可靠的海底地形統計信息。

3 結論

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

2)配置法的關鍵是計算轉換函數，這需要計算海底地形與重力異常之間的相關函數，因而需要預先已知海底地形較準確的統計信息，這與海底地形反演相矛盾，實際中較少采用。

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

4)垂直重力梯度異常數據也可用于聯合船測海深數據計算海底地形模型，且計算精度與重力異常反演一致。

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，Forsberg 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 黃謨濤.利用衛星測高資料反演海底地形研究［J］.武漢大學學報:信息科學版，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 吳云孫.利用測高重力梯度異常反演中國南海海底地形［J］.武漢大學學報:信息科學版，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 李大煒.利用衛星測高數據反演中國海及鄰近海域海底地形［J］.大地測量與地球動力學，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］.大地測量與地球動力學，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，Frequency 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

中國地震局地震研究所所長基金重點項目(IS201326125);國家測繪地理信息局測繪基礎研究基金項目(13－01－01);國家自然科學基金項目(41204019，41304003)。

胡敏章，男，1985年生，博士，助理研究員，研究方向為海底地形與重力均衡。E－mail:huminzhang@126.com。