基于衛(wèi)星測(cè)高數(shù)據(jù)的海底地形反演及分析

眭曉虹 張潤寧 萬曉云 李洋

(1 航天東方紅衛(wèi)星有限公司，北京 100094) (2 錢學(xué)森空間技術(shù)實(shí)驗(yàn)室，北京 100094) (3 北京空間飛行器總體設(shè)計(jì)部，北京 100094)

?

基于衛(wèi)星測(cè)高數(shù)據(jù)的海底地形反演及分析

眭曉虹1,2張潤寧3萬曉云2李洋1

(1 航天東方紅衛(wèi)星有限公司，北京 100094) (2 錢學(xué)森空間技術(shù)實(shí)驗(yàn)室，北京 100094) (3 北京空間飛行器總體設(shè)計(jì)部，北京 100094)

由于受限于觀測(cè)條件，海底地形模型在分辨率與精度兩方面較陸地和大地水準(zhǔn)面模型有很大的差距。傳統(tǒng)的海底地形測(cè)量是利用聲吶船測(cè)技術(shù)獲取水深數(shù)據(jù)，但該測(cè)量方法測(cè)量效率低、測(cè)量點(diǎn)分布稀疏不均勻。文章詳細(xì)研究并推導(dǎo)了基于衛(wèi)星測(cè)高數(shù)據(jù)進(jìn)行海底地形反演的方法，采用衛(wèi)星測(cè)高重力異常數(shù)據(jù)和少量船測(cè)水深數(shù)據(jù)，對(duì)北太平洋部分區(qū)域進(jìn)行了海底地形反演，得到1′×1′分辨率的海底地形模型，并對(duì)反演關(guān)系中的非線性項(xiàng)進(jìn)行了分析。仿真結(jié)果驗(yàn)證了此方法的可行性，指出海底地形與重力異常數(shù)據(jù)在一定波段內(nèi)以線性關(guān)系為主，疊加非線性二次項(xiàng)和三次項(xiàng)可較好地提升反演精度，高于三次項(xiàng)的部分對(duì)反演的結(jié)果沒有明顯的作用，可以忽略。

海底地形；衛(wèi)星測(cè)高；重力異常；聲吶船測(cè)技術(shù)

1 引言

海底地形是描述海洋水下固體表面起伏的物理概念，對(duì)海底地形的精確探測(cè)不僅對(duì)研究地球物理科學(xué)具有重大的意義，同時(shí)對(duì)潛艇水下導(dǎo)航等軍事活動(dòng)、海底管線鋪設(shè)等經(jīng)濟(jì)活動(dòng)都具有指導(dǎo)作用[1-2]。在傳統(tǒng)的海洋測(cè)繪中，聲吶船測(cè)技術(shù)是獲取測(cè)量船到海底地形起伏面的深度(即水深)數(shù)據(jù)的核心手段。利用測(cè)量船搭載聲吶技術(shù)測(cè)量水深具有精度高的優(yōu)點(diǎn)，但缺點(diǎn)是效率低、測(cè)量點(diǎn)分布稀疏不均勻，無法有效地實(shí)現(xiàn)高分辨率、高精度、全球化海底地形的探測(cè)。

自20世紀(jì)70年代發(fā)射第一顆搭載雷達(dá)高度計(jì)的測(cè)高衛(wèi)星起，應(yīng)用衛(wèi)星測(cè)高技術(shù)獲取全球海洋大地水準(zhǔn)面高的數(shù)據(jù)變?yōu)榱丝赡躘3-5]。大地水準(zhǔn)面作為最接近全球平均海面的重力場(chǎng)等勢(shì)面，其形狀反映了地球內(nèi)部物質(zhì)構(gòu)造、密度及分布等信息，與海底地形在一定程度上具有高相關(guān)性[6]。Dixon、Ramillien、Hwang、黃謨濤等多名國內(nèi)外研究人員基于衛(wèi)星測(cè)高數(shù)據(jù)提取海面重力測(cè)量數(shù)據(jù)，并對(duì)部分海域進(jìn)了海底地形反演[7-10]，但這些反演過程只考慮海面重力數(shù)據(jù)與海底地形線性相關(guān)性，對(duì)于非線性的相關(guān)程度未予考慮。

本文基于Smith&Sandwell(S&S)的經(jīng)典方法[1]，系統(tǒng)探討了應(yīng)用衛(wèi)星測(cè)高數(shù)據(jù)進(jìn)行海底地形反演的方法，并進(jìn)行了仿真分析。仿真結(jié)果表明，基于S&S方法可以有效地反演海底地形。針對(duì)前人忽略的反演關(guān)系中非線性項(xiàng)部分，文中亦做了仿真分析，結(jié)果表明，海底地形與重力異常數(shù)據(jù)在一定波段內(nèi)以線性關(guān)系為主，疊加非線性二次項(xiàng)和三次項(xiàng)可較好地提升反演精度，高于三次項(xiàng)的部分對(duì)反演的結(jié)果沒有明顯的作用。

2 衛(wèi)星測(cè)高數(shù)據(jù)反演海底地形的原理及方法

2.1 衛(wèi)星測(cè)高數(shù)據(jù)計(jì)算海洋重力異常

基于Sandwell方法，可從衛(wèi)星測(cè)高數(shù)據(jù)得到沿軌海面高程斜率，結(jié)合測(cè)高點(diǎn)的位置和時(shí)間信息得到沿軌大地水準(zhǔn)面起伏對(duì)空間的導(dǎo)數(shù)，從而計(jì)算海面垂線偏差。如圖1所示，以相對(duì)于平均海面z0深度的參考水深面作xy平面建立直角坐標(biāo)系，x軸指向地理東向，y軸指向地地理北向，z軸與x和y軸構(gòu)成右手坐標(biāo)系。記海面東向和北向垂線偏差的分量分別為

(1)

(2)

式中：N(x,y)表示點(diǎn)(x,y)處大地水準(zhǔn)面高度，即大地水準(zhǔn)面與參考橢球面的距離差。

圖1 平均海面與海底地形坐標(biāo)示意圖Fig.1 Coordinate system for mean sea surface and sea floor topography

根據(jù)逆Vening-Meinesz公式，可由垂線偏差計(jì)算海洋重力異常Δg[11]為

η(P)sinα]dσ

(3)

式中：g0為海平面的重力加速度，通常取9.8m·s-2；ξ(P)和η(P)分別為P處垂線偏差的北向和東向分量；α為流動(dòng)點(diǎn)至測(cè)量點(diǎn)的方位角；H′(ψ)為積分核函數(shù)。

(4)

式中：ψ為流動(dòng)點(diǎn)到計(jì)算點(diǎn)的球面角距。

至此，已建立測(cè)高數(shù)據(jù)與海洋重力異常的相互關(guān)系。根據(jù)物理大地測(cè)量學(xué)相關(guān)理論，海底地形起伏變化將引起局部重力場(chǎng)擾動(dòng)，記為擾動(dòng)重力δg，該重力場(chǎng)擾動(dòng)包含于海洋重力異常。針對(duì)重力異常數(shù)據(jù)在相關(guān)波段內(nèi)濾波，保留與海底地形具有高相關(guān)性的擾動(dòng)重力，才可進(jìn)而反演海底地形。

2.2 利用擾動(dòng)重力反演海底地形

假設(shè)存在z=h(r)的一個(gè)海底地形起伏面，如圖1所示，h(r)為相對(duì)于參考水深面的海底地形高度。圖中海面上P點(diǎn)處海底地形產(chǎn)生的擾動(dòng)位記為

(5)

式中：G=6.672 6×10-11Nt·m2·kg-2為地球引力常數(shù)；r0和r分別為計(jì)算點(diǎn)的位置矢量和海底地形起伏面z=h(r)的流動(dòng)體元的位置矢量；Δρ為海底固體巖石層與海水的密度差。

采用頻域分析法，較之經(jīng)典的時(shí)域分析法可以快速地體現(xiàn)系統(tǒng)的響應(yīng)過程，高效率地計(jì)算得到系統(tǒng)響應(yīng)輸出，并為相關(guān)波段的截取提供便利。將式(5)在z=z0平面內(nèi)進(jìn)行二維傅里葉變換可快速得到

(6)

海面擾動(dòng)重力是重力擾動(dòng)位U(r0)在垂直方向上的導(dǎo)數(shù)，計(jì)算式(6)并根據(jù)傅里葉變換的性質(zhì)可得頻域內(nèi)擾動(dòng)重力δg的表達(dá)式：

F[δg]=|k|F[U(r0)]=

(7)

根據(jù)式(7)即可利用擾動(dòng)重力δg數(shù)據(jù)反演海底地形高度h(r)[1,12-13]。

2.3 反演方法

本文基于式(7)，利用規(guī)則的衛(wèi)星測(cè)高重力異常數(shù)據(jù)結(jié)合稀疏的船測(cè)水深數(shù)據(jù)反演海底地形。對(duì)測(cè)高重力異常數(shù)據(jù)在與海底地形有高相關(guān)性的波段進(jìn)行濾波，得到擾動(dòng)重力。文中以部分船測(cè)點(diǎn)作為控制點(diǎn)，將控制點(diǎn)處的水深數(shù)據(jù)作為真實(shí)水深數(shù)據(jù)，濾除與擾動(dòng)重力相關(guān)性差的長波水深數(shù)據(jù)得到殘余水深。以船測(cè)殘余水深數(shù)據(jù)與船測(cè)航跡上的擾動(dòng)重力建立相關(guān)關(guān)系，得到只考慮一階項(xiàng)的線性系數(shù)和考慮多階項(xiàng)的非線性系數(shù)。在船測(cè)航跡之外的區(qū)域，將得到的線性和非線性系數(shù)疊加到擾動(dòng)重力的作用中，從而得到非船測(cè)航跡處預(yù)測(cè)水深。結(jié)合預(yù)測(cè)水深數(shù)據(jù)與長波水深數(shù)據(jù)，即可得到計(jì)算區(qū)域內(nèi)的全波段水深，即海底地形。將未參與計(jì)算的船測(cè)點(diǎn)作為外部校核點(diǎn)，用以衡量海底地形的反演精度。

3 數(shù)據(jù)準(zhǔn)備

仿真實(shí)驗(yàn)針對(duì)區(qū)域40°N～50°N，150°E～160°E的北太平洋海域內(nèi)，基于一組衛(wèi)星測(cè)高重力異常數(shù)據(jù)進(jìn)行海底地形反演，并應(yīng)用一組稀疏船測(cè)水深數(shù)據(jù)作為控制點(diǎn)輔助仿真。重力異常數(shù)據(jù)來源于丹麥科技大學(xué)空間研究所，數(shù)據(jù)分辨為1′×1′，數(shù)據(jù)總點(diǎn)數(shù)為361 201，該數(shù)據(jù)由多年的測(cè)高衛(wèi)星數(shù)據(jù)計(jì)算得到，數(shù)據(jù)如圖2所示。

圖2 重力異常數(shù)據(jù)Fig.2 Gravity anomalies data

稀疏船測(cè)數(shù)據(jù)來源于美國大氣海洋管理局(NOAA)，船測(cè)數(shù)據(jù)總點(diǎn)數(shù)為57 986。稀疏船測(cè)水深數(shù)據(jù)(如圖3所示)表示測(cè)量船(海面)到海底地形起伏面的深度，取向上為正。

圖3 船測(cè)水深數(shù)據(jù)(圖中包含測(cè)量船的航跡線)Fig.3 Shipboard bathymetry data(with survey ship routes)

文中取56 826個(gè)船測(cè)數(shù)據(jù)點(diǎn)作為控制點(diǎn)帶入海底地形反演的計(jì)算，數(shù)據(jù)點(diǎn)分布如圖4所示。

圖4 船測(cè)水深控制點(diǎn)數(shù)據(jù)(無船測(cè)航跡線)Fig.4 Shipboard bathymetry data of reference points(without survey ship routes)

其余1160個(gè)船測(cè)數(shù)據(jù)點(diǎn)作為核驗(yàn)點(diǎn)帶入海底地形反演精度的核驗(yàn)，數(shù)據(jù)點(diǎn)分布如圖5所示。

圖5 船測(cè)水深核驗(yàn)點(diǎn)數(shù)據(jù)Fig.5 Shipboard bathymetry data of check points

針對(duì)船測(cè)水深數(shù)據(jù)中部分?jǐn)?shù)據(jù)偏差較大的情況，需要對(duì)船測(cè)水深數(shù)據(jù)進(jìn)行粗差校正。為了方便粗差校正，并濾除長波水深數(shù)據(jù)，首先對(duì)不規(guī)則的船測(cè)水深數(shù)據(jù)插值得到格網(wǎng)化的原始船測(cè)水深，如圖6所示；經(jīng)過粗差校正后的水深數(shù)據(jù)如圖7所示，從圖6與圖7的對(duì)比可看出，粗差校正有效地剔除船測(cè)水深中偏差較大的數(shù)據(jù)，從而提高仿真實(shí)驗(yàn)的數(shù)據(jù)可信性。

圖6 船測(cè)水深格網(wǎng)化后的三維圖像Fig.6 3D presentation for gridded shipboard bathymetry

圖7 粗差校正后的船測(cè)水深三維圖像Fig.7 3D presentation for gridded shipboard bathymetry after cross error correction

4 仿真分析

仿真基于文中第2節(jié)的反演原理和方法，利用第3節(jié)的實(shí)驗(yàn)準(zhǔn)備數(shù)據(jù)進(jìn)行仿真實(shí)驗(yàn)，具體的仿真流程如圖8所示。

圖8 仿真實(shí)驗(yàn)流程圖Fig.8 Flow chart for simulation

根據(jù)重力異常與水深在海底地形波長為20～200 km波段內(nèi)具有高相關(guān)性的理論，對(duì)船測(cè)水深數(shù)據(jù)進(jìn)行低通濾波，濾波器設(shè)計(jì)如圖9所示，可獲得波長大于200 km的長波船測(cè)水深數(shù)據(jù)，如圖10所示(圖中將船測(cè)水深轉(zhuǎn)化為地形高度)；同時(shí)對(duì)重力異常數(shù)據(jù)進(jìn)行帶通濾波，濾波器設(shè)計(jì)如圖11所示，可獲得波長介于20 km與200 km之間的擾動(dòng)重力數(shù)據(jù)并向下延拓至參考水深處(文中參考水深為z0=4 km)，如圖12所示。

圖9 低通濾波器Fig.9 Low pass filter

圖10 長波船測(cè)海底地形數(shù)據(jù)Fig.10 Long wavelength shipboard bathymetry

圖11 帶通濾波器Fig.11 Band pass filter

針對(duì)船測(cè)殘余水深數(shù)據(jù)與20～200 km的船測(cè)航跡上的擾動(dòng)重力數(shù)據(jù)進(jìn)行回歸分析并均勻插值，得到船測(cè)航跡上海底地形與擾動(dòng)重力的線性系數(shù)和多階非線性系數(shù)。將仿真區(qū)域的擾動(dòng)重力與相關(guān)系數(shù)相乘即可得到短波水深數(shù)據(jù)，將計(jì)算水深與長波水深相加即可得到全區(qū)域的反演水深。

圖12 帶通濾波并向下延拓至參考水深處的擾動(dòng)重力數(shù)據(jù)Fig.12 Gravity anomaly disturbance data being band pass filtered and downward continued to reference water depth

經(jīng)過計(jì)算，最終得到1′×1′的海底地形模型。圖13所示為僅考慮一階線性項(xiàng)的海底地形反演結(jié)果，地形高度以相對(duì)于參考水深面(低于平均海面4 km處的參考面)的地形起伏高度進(jìn)行表示。由于考慮非線性項(xiàng)與只考慮線性項(xiàng)的反演結(jié)果差異遠(yuǎn)小于海底地形的量級(jí)，因此考慮非線性項(xiàng)與只考慮線性項(xiàng)的反演結(jié)果圖沒有明顯區(qū)別，文中只給出取線性項(xiàng)的海底地形反演圖，具體非線性項(xiàng)的影響在表1中以數(shù)據(jù)形式體現(xiàn)。

圖13 取線性項(xiàng)的海底地形反演結(jié)果Fig.13 Predicted bathymetry with first order

反演的精度可由反演水深結(jié)果與未參與計(jì)算的船測(cè)控制點(diǎn)水深進(jìn)行比對(duì)，表1對(duì)考慮一階線性項(xiàng)和多階非線性項(xiàng)的反演數(shù)據(jù)精度進(jìn)行統(tǒng)計(jì)。

表1 海底地形反演數(shù)據(jù)統(tǒng)計(jì)

將反演得到的海底地形與船測(cè)數(shù)據(jù)進(jìn)行相關(guān)性檢驗(yàn)可得數(shù)據(jù)如圖14所示。

圖14 反演海底地形與船測(cè)數(shù)據(jù)的相關(guān)性Fig.14 Coherence between predicted bathymetry and shipboard bathymetry

根據(jù)表1及圖14可看出，在20～200 km的波段范圍內(nèi)，重力異常與海底地形以線性關(guān)系為主，疊加非線性二次項(xiàng)和三次項(xiàng)將反演精度的平均值量級(jí)提升至米級(jí)，反演得到的海底地形與船測(cè)數(shù)據(jù)的相關(guān)性也有較大提高，高于三次項(xiàng)的部分對(duì)反演的結(jié)果沒有明顯的作用，同時(shí)由于控制點(diǎn)數(shù)據(jù)量有限，更高階的關(guān)系系數(shù)并不能準(zhǔn)確地獲得，從而導(dǎo)致從五階之后階數(shù)增大反演結(jié)果卻并不理想。

5 結(jié)束語

本文詳細(xì)研究并推導(dǎo)了由衛(wèi)星測(cè)高數(shù)據(jù)反演海底地形的原理及方法，基于此方法，利用衛(wèi)星測(cè)高重力異常數(shù)據(jù)和船測(cè)水深數(shù)據(jù)對(duì)北太平洋10°×10°的區(qū)域進(jìn)行海底地形反演。數(shù)據(jù)及分析結(jié)果表明：該方法可以有效地反演海底地形，且海底地形與重力異常數(shù)據(jù)在20～200 km波段內(nèi)主要呈現(xiàn)線性關(guān)系，非線性二次項(xiàng)和三次項(xiàng)對(duì)反演精度的提升具有一定的作用，高于三次項(xiàng)的部分對(duì)反演的結(jié)果沒有明顯的作用。由于波長小于20 km和波長大于200 km海底地形起伏引起重力異常的物理機(jī)理非常復(fù)雜，且重力異常在向下延拓時(shí)會(huì)出現(xiàn)解不唯一且不穩(wěn)定的問題，由此產(chǎn)生的低于20 km和高于200 km波段的重力異常數(shù)據(jù)不能準(zhǔn)確地與船測(cè)水深數(shù)據(jù)建立關(guān)系，也是導(dǎo)致最終海底地形反演精度不高的原因之一。因此，未來對(duì)全波段重力異常與水深關(guān)系的研究將具有非常重要的意義。

References)

[1]Walter H F Smith, David T Sandwell. Bathymetry prediction from dense satellite altimetry and sparse shipboard bathymetry[J]. Journal of Geophysical Research, 1994, 99(B11): 21803-21824

[2]T H Bell. Statistical features of sea-floor topography[J]. Deep Sea Research, 1975, 22: 883-892

[3]D C McAdoo, K M Marks. Gravity fields from Geosat and ERS-1 satellite altimetry for geodynamic applications[J]. Eos Trans, 1992(74): 133

[4]D T Sandwell, Walter H F Smith. Global marine gravity from ERS-1, Geosat and Seasat reveals new tectonic fabric[J]. Eos Trans,1992, AGU(73):133

[5]G A Neumann, D W Forsyth, D T Sandwell. Comparison of marine gravity from shipboard and high-density satellite altimetry along the mid-Atlantic ridge, 30.5°-35.5°S[J]. Geophysical Research Letters, 1993, 20(15): 1639-1642

[6]呂華慶. 物理海洋學(xué)基礎(chǔ)[M]. 北京：海洋出版社, 2012

LU Huaqing. Fundamental of physical oceanography[M]. Beijing: Ocean Press, 2012 (in Chinese)

[7]T H Dixon. Bathymetric predition from Seasat altimeter data[J]. Journal of Geophysics and Research, 1983, 88: 1563-1571

[8]G Ramillien, A cazenave. global bathymetry derived from altimeter data of the ERS-1 geodetic mission[J]. Journal of Geodynamics, 1997, 23: 129-149

[9]C Hwang. A bathymetric model for the South China Sea from satellite altimetry and depth data[J]. Marine Geodesy, 1999, 22(1): 37-51

[10]黃謨濤,翟國君，歐陽永忠, 等. 衛(wèi)星測(cè)高資料在反演海底地形中的應(yīng)用[J]. 海洋測(cè)繪, 2002, 22(1): 3-7

Huang Motao, Zhai Guojun, Ouyang Yongzhong, et al. Applications of satellite altimetry data in bathymetry prediction[J]. Hydrographic Surveying and Charting, 2002, 22(1): 3-7 (in Chinese)

[11]C Hwang. Inverse Vening-Meinesz formula and deflection-geoid formula: applications to the predictions of gravity and geoid over the South China Sea[J]. Journal of Geodesy, 1998, 72(5): 304-312

[12]Hu Minzhang, Li Jiancheng, Li Hui, et al. Bathymetry predicted from vertical gravity gradient anomalies and ship soundings[J]. Geodesy and Geodynamics, 2014, 5(1): 41-46

[13]歐陽明達(dá),孫中苗,翟振和, 等.采用重力異常的導(dǎo)納理論推估海底地形[J]. 測(cè)繪學(xué)報(bào), 2015, 44(10): 1092-1099

Ouyang Mingda, Sun Zhongmiao, Zhai Zhenhe, et al. Bathymetry prediction based on the admittance theory of gravity anomalies[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(10): 1092-1099 (in Chinese)

(編輯：張小琳)

Bathymetry Prediction and Analysis Based on Satellite Altimetry Data

SUI Xiaohong1,2ZHANG Running3WAN Xiaoyun2LI Yang1

(1 DFH Satellite Co., Ltd., Beijing 100094, China) (2 Qian Xuesen Laboratory of Space Technology, Beijing 100094, China) (3 Beijing Institute of Spacecraft System Engineering, Beijing 100094, China)

Both of the resolution and accuracy of bathymetry model is far from the model of the earth land and geoid, because of observing difficulty. The traditional measurement of the bathymetry is using shipboard sounding technique, which has the limitations of low efficiency, sparse and nonuniform covering. In this paper, bathymetry prediction method based on satellite altimetry data is presented. An 1′×1′ bathymetry model of a region in north Pacific is predicted through simulation in this paper, and the nonlinear relationship between satellite altimetry data and bathymetry data is analyzed. The results show that, it is feasible to use this method predicting bathymetry. The relationship between satellite altimetry gravity anomalies and bathymetry is mainly linear in certain wavelength; adding the second and third order nonlinear terms will improve the accuracy of the inversion, however the higher order terms have little effect.

bathymetry; satellite altimetry; gravity anomalies; shipboard sounding technique

2017-02-01；

2017-05-26

國家自然科學(xué)基金(41404019，41674026)，中國科學(xué)院太空應(yīng)用重點(diǎn)實(shí)驗(yàn)室開放基金(CSU-WX-A-KJ-2016-044)

眭曉虹，女，博士研究生，研究方向?yàn)樾l(wèi)星測(cè)高總體設(shè)計(jì)及應(yīng)用。Email：yuanyuansui@hotmail.com。

P228

A

10.3969/j.issn.1673-8748.2017.03.020