陽凡林，沈瑞杰，梅賽，屠澤杰，辛明真

(1.山東科技大學 測繪與空間信息學院，山東 青島 266590；2.自然資源部海洋測繪重點實驗室，山東 青島 266590；3.中國地質調查局 青島海洋地質研究所，山東 青島 266237)

1 引言

海洋約占地球表面積的71%，蘊藏著極其豐富的礦產和生物資源。海底地形數(shù)據(jù)是研究海洋生物、海洋化學及海洋地球物理學等相關學科的基礎。自20 世紀50 年代回波測深技術廣泛應用以來，搭載單波束和多波束的研究船可以開展高精度、高分辨率的海底地形勘測。然而船載水深測量存在效率低、測量結果分布不均勻的問題，使用回聲測深技術對全球海底地形進行完全覆蓋式勘測仍需要很長時間[1]。

衛(wèi)星測高技術的出現(xiàn)，可以全天候、連續(xù)準確地在全球范圍內直接獲取海面高度的觀測數(shù)據(jù)，其測量效率和覆蓋范圍是船載水深測量無法比擬的。隨著衛(wèi)星測高資料的積累，利用測高數(shù)據(jù)獲取的海洋重力信息已成為反演海底地形的主要數(shù)據(jù)來源。目前，利用重力數(shù)據(jù)進行海底地形反演的研究方法主要有：重力地質法（Gravity Geologic Method，GGM）[2–3]、導納函數(shù)法（Admittance Method）[4]、最小二乘配置法（Least Squares Collocation，LSC）[5–6]和Smith and Sandwell（SAS）法[7]。其中，GGM 和SAS 法使用最為廣泛，兩者大多以重力異常數(shù)據(jù)作為輸入。Kim 等[8]將GGM 與SAS法構建的海底地形模型進行比較，發(fā)現(xiàn)GGM 在預測短波長（≤12 km）分量方面具有優(yōu)勢，而SAS 方法在預測長波長（≥25 km）分量方面具有優(yōu)勢；李倩倩和鮑李峰[9]，彭聰?shù)萚10]均對SAS 方法和GGM 進行了比較分析，后者給出了兩種方法的水深和海底地形復雜程度的適用條件；范雕等[11]應用導納函數(shù)法建立了不同均衡補償模式下的海底地形模型，并對各項地球物理參數(shù)的影響進行了詳細分析；郭金運等[12]利用迭代延拓法，提升了重力異常向下延拓的精度。

除重力異常外，重力垂直梯度異常（Vertical Gravity Gradient Anomalies，VGG）也是反演海底地形的重要數(shù)據(jù)，其反映了重力場在垂直方向上的變化情況，且重力垂直梯度異常相比于重力異常含有更多的高頻信號成分，能夠更好地描述異常特征[13–14]；Wang[15]指出大洋地殼的質量補償對重力梯度的影響很小，提出可以利用重力垂直梯度異常反演水深的思路；Wessel 和Lyons[16]通過對高斯型海山模型的研究，指出重力垂直梯度異常具有放大短波信號、抑制長波信號的作用；吳云孫等[17]推導了重力垂直梯度異常與海底地形之間的關系并反演了南海海域的海底地形模型；胡敏章等[18]通過疊加重力異常（波長<100 km）和重力垂直梯度異常（波長100～200 km）反演的海底地形構建了中國海域的海底地形模型，其精度與V15.1 相近；范雕等[19]利用多元回歸分析技術，通過解算重力異常和重力垂直梯度異常的偏回歸系數(shù)和常數(shù)項來反演海底地形。

目前，聯(lián)合重力異常和重力垂直梯度異常數(shù)據(jù)進行海底地形反演的研究較少。本文針對重力異常、重力垂直梯度異常和船測水深數(shù)據(jù)各自的特性，通過穩(wěn)健回歸分析降低了船測數(shù)據(jù)和重力數(shù)據(jù)中異常點的偏差影響，分別使用重力異常和重力垂直梯度異常數(shù)據(jù)反演了太平洋皇帝山海域（35°～45°N，165°～175°E）的海底地形模型，采用滑動窗口賦權的方法，選擇滑動窗口內的最優(yōu)權值，聯(lián)合兩種重力數(shù)據(jù)構建的模型。最后，以船測水深數(shù)據(jù)為外部檢核，與DTU18 和SIO V23.1 全球海底地形模型以及地形剖面處的實測水深數(shù)據(jù)進行比較分析。

2 研究區(qū)域和數(shù)據(jù)來源

夏威夷皇帝海山鏈位于太平洋西北側，該區(qū)域存在海底山脈且兩側海隆聚集，海底深度跨度大（200～8 000 m），是開展海底地形反演的熱點區(qū)域。本文選取的研究區(qū)域位于其南部（35°～45°N，165°～175°E），考慮到邊緣效應，將研究區(qū)域的經緯度均向外擴充1°，研究區(qū)域如圖1 所示。

圖1 研究區(qū)域與船載測深軌跡分布Fig.1 Study area and shipborne bathymetry trace distribution map

2.1 船載測深數(shù)據(jù)

本文采用的船載測深數(shù)據(jù)來自于美國國家地球物理數(shù)據(jù)中心（National Geophysical Data Center，NGDC）發(fā)布的實測數(shù)據(jù)，其中單波束數(shù)據(jù)時間跨度為42 a（1961–2002 年），多波束數(shù)據(jù)時間跨度為16 a（2004–2019 年）。由于早期數(shù)據(jù)可能存在測量誤差，以全球地形模型ETOPO1 作為約束條件進行互差比較，依照“3σ”準則，采用10′×10′的窗口大小，5′的移動步長滑動剔除單波束數(shù)據(jù)中的粗差點，得到61 169 個單波束船測點。均勻選取其中4/5 的船測點作為控制點，剩余船測單波束數(shù)據(jù)和濾波抽稀后的多波束數(shù)據(jù)用于后續(xù)模型的精度評定及地形剖面對比。

2.2 重力數(shù)據(jù)

重力異常數(shù)據(jù)（圖2a）和重力垂直梯度異常數(shù)據(jù)（圖2b）均來源于加州大學圣地亞哥分校（University of California,San Diego，UCSD）的斯克里普斯海洋研究所（Scripps Institution of Oceanography，SIO）于2021年8 月發(fā)布的V31.1 版本，格網分辨率為1′×1′。對比先前版本，V31.1 模型補充了新的Altika、Cryosat LRM、Cryosat SAR 和Sentinel-3A/B 衛(wèi)星測高數(shù)據(jù)。為了方便說明，下文將重力異常和重力垂直梯度異常數(shù)據(jù)合稱為重力數(shù)據(jù)。

圖2 重力模型示意圖Fig.2 Gravimetric model diagram

2.3 海底地形模型

本文采用目前廣泛使用的ETOPO1、DTU18 和SIO V23.1 海底地形模型（1′×1′）作為參考模型。

ETOPO1 由美國國家地球物理中心（National Geophysical Data Center，NGDC）和美國國家海洋和大氣管理局（National Oceanic and Atmospheric Administration，NOAA）于2008 年8 月發(fā)布，該模型集成了來自區(qū)域和全球數(shù)據(jù)集的地形、測深和海岸線數(shù)據(jù)。

DTU18 由丹麥科技大學（Technical University of Denmark，DTU）于2019 年9 月發(fā)布，該模型補充使用了3 a 的Sentinel-3A 衛(wèi)星和7 a 的Cryosat-2 LRM 衛(wèi)星數(shù)據(jù)，并使用FES2014 海潮模型進行潮汐改正。

SIO V23.1 由美國SIO 的Sandwell 團隊于2021 年11 月發(fā)布，該模型結合了最新的重力模型數(shù)據(jù)、衛(wèi)星測高數(shù)據(jù)以及優(yōu)化的船測水深數(shù)據(jù)，是公認精度最高的全球海底地形模型之一。

3 海底地形反演方法

Smith 和Sandwell[7]于1994 年提出了依據(jù)線性回歸分析技術反演海底地形的方法。其基本原理是利用重力數(shù)據(jù)和海底地形在強相關性波段內的線性關系，反演一定波段的海底地形，并疊加其他波段的船測水深數(shù)據(jù)得到全波段海底地形模型。由于線性回歸分析技術會受到重力數(shù)據(jù)或船測數(shù)據(jù)中一定異常點的影響，因此采用穩(wěn)健回歸分析法，提高線性回歸分析中比例因子的魯棒性；并提出滑動窗口賦權，聯(lián)合由重力異常和重力垂直梯度異常構建的模型，通過對窗口內的兩種模型分配最優(yōu)權值，構建研究海域最終的海底地形模型。下面對使用的方法進行具體介紹。

3.1 反演波段的選擇

Parker[20]于1973 年將快速傅里葉變換[21]應用到位理論，根據(jù)Parker 公式，在不顧及地殼撓曲均衡影響下，海底地形與重力異常G0(k)之間的關系可以描述為

同理，海底地形與重力垂直梯度異常G1(k)之間的關系可以描述為[18]

式中，G為地球引力常數(shù)；k為徑向頻率(kx,ky)=(1/λx,1/λy) ），其中 (kx,ky)和(λx,λy)分別為x和y方向的頻率和波長；d為平均海深值；ρc和ρw分別為地殼和海水密度；Z(k)為導納函數(shù)。

由Parker 公式可知，在頻率域內重力數(shù)據(jù)和海底地形在一定波段范圍內具有相干性，若將重力數(shù)據(jù)和水深看作兩種不同的信號，兩者的相干性可以表示為

式中，r2(k)為互譜相干函數(shù)，當r2(k)趨近于1 時，說明兩種信號呈線性相關，趨近于0 時，說明兩種信號完全不相干；G(k)、H(k)分別代表重力信號和地形信號的傅里葉變換；G*(k)、H*(k)表示為G(k)、H(k)的復數(shù)共軛。研究海域內重力數(shù)據(jù)與海底地形相干性關系如圖3 所示。

圖3 重力數(shù)據(jù)與海底地形的相干性Fig.3 Coherence between gravimetric data and seafloor topography

重力垂直梯度異常在50～200 km 波長范圍內與海底地形相干性大于0.5，而重力異常數(shù)據(jù)在50～400 km 波長范圍內相干性大于0.5，考慮到在波長大于200 km 時，會受到地殼均衡補償影響[22]，海底地形產生很少或沒有重力場變化，因此本文選擇的反演波段為50～200 km。

3.2 穩(wěn)健回歸分析

穩(wěn)健回歸分析的原理是基于最小二乘估計的迭代加權法，通過對權的優(yōu)選以提高擬合精度與可靠性，使參數(shù)的估值盡可能避免粗差的影響，得到正常模式下的最佳估值[23]，對于線性回歸模型

Huber 函數(shù)將殘差分為兩組：一組殘差超過初始標準差的某個因子c，另一組殘差不超過這個限制。u代表標準化的殘差指標（ui=vi/σ0），σ0為單位權中誤差，c表示調和系數(shù)，本文取c=2。參數(shù)a的穩(wěn)健估計可用迭代加權最小二乘求解，陳艷國[25]指出了具體求解步驟。反演波段的重力數(shù)據(jù)和船測水深數(shù)據(jù)線性擬合結果如圖4 所示。

圖4 反演波段重力數(shù)據(jù)與殘余海深線性擬合結果Fig.4 Results of linear fitting between gravity data of inversion band and residual sea depth

可以看出，在研究海域，重力數(shù)據(jù)與殘余海深在反演波段表現(xiàn)出了良好的線性關系。文獻[26]指出，傳統(tǒng)的最小二乘分析相比于穩(wěn)健回歸分析更容易受到離群點的影響，導致其斜率（即比例因子）傾向于離群點，而利用Huber 損失函數(shù)的穩(wěn)健回歸分析通過給予離群點較小的權值，能夠減小對離群點的敏感度問題，以此提高比例因子的魯棒性。

3.3 滑動窗口賦權

滑動窗口賦權是建立在最小二乘賦權[27]的基礎上，將重力異常反演的模型G0(i)和重力垂直梯度異常反演的模型G1(i)與船測水深值H(i)進行計算。假設滑動窗口內G0(i)所占的權重為xi，則G1(i)所占的權重為1?xi，融合得到的G′(i)可表示為

將融合后的G′(i)和H(i)進行比較，以其差值的標準差最小為準則分配最優(yōu)權值，計算公式為

以此對滑動窗口內的G0(i)和G1(i)分配最優(yōu)權值，進而融合得到區(qū)域不同地形起伏條件下的權值分配。研究海域內，重力垂直梯度異常模型的權值分配如圖5 所示。

圖5 重力垂直梯度異常模型權值xi 分配示意圖，相應的，重力異常模型的權值分配為1?xiFig.5 The schematic diagram of weight distribution xi of vertical gravity gradient anomalies model.Accordingly,the weight distribution of gravity anomalies model is 1?xi

基于上述分析，本文反演海底地形模型的基本流程如圖6 所示。

圖6 海底地形反演流程Fig.6 Flowchart of seafloor topography inversion

①通過相干性分析，獲取重力數(shù)據(jù)和海底地形之間相干性高的波段作為反演波段，將重力數(shù)據(jù)進行帶通濾波和向下延拓處理，獲得反演波段的重力數(shù)據(jù)G(x)并插值到船測控制點處；

② 將格網化的船載測深數(shù)據(jù)進行低通濾波得到長波海深模型，內插至船測控制點處與海深值做差得到控制點處殘余海深；

③將步驟①的中波重力數(shù)據(jù)與步驟②的殘余海深以滑動窗口方式進行穩(wěn)健回歸分析，解算格網點處比例因子S(x)和常數(shù)項；

④ 將比例因子S(x)和重力數(shù)據(jù)G(x)相乘并疊加常數(shù)項得到反演波段海底地形，與長波海深模型相加得到中長波海深模型；

⑤ 將中長波海深模型內插至船測控制點得到控制點處中長波水深值，與控制點處水深測量值做差得到控制點處短波水深并格網化，最后與中長波海底地形模型相加得到全波段海底地形模型；

⑥ 采用滑動窗口賦權融合重力異常和重力垂直梯度異常數(shù)據(jù)構建的海底地形模型，得到最終的海底地形模型。

4 模型構建結果及精度評價

4.1 多源海底地形模型的構建

以20′×20′的窗口大小[28]，通過穩(wěn)健回歸分析擬合格網處比例因子和常數(shù)項。分別將重力異常和重力垂直梯度異常構建的模型稱為SGA（Single Gravity Anomalies）和SVG（Single Vertical Gradient）。以SGA和SVG 模型與船測檢核點間標準差最小為準則，采用滑動窗口賦權法以30′×30′的窗口大小，將SGA 和SVG 模型分配窗口內的最優(yōu)權值，融合構建的最終模型稱為MGM（Multiple Gravity Model），反演結果如圖7 所示。

圖7 MGM 海底地形模型Fig.7 The seafloor topography model of multiple gravity model

4.2 與單一重力數(shù)據(jù)源模型比較分析

為驗證本文方法的優(yōu)越性，對SGA、SVG、傳統(tǒng)最小二乘賦權法（MGM_LS）和MGM 4 種模型進行比較分析，考慮到構建以上4 種模型時均未使用多波束數(shù)據(jù)進行約束，且多波束數(shù)據(jù)集中分布在地形起伏變化劇烈海域，因此使用該數(shù)據(jù)作為外部檢核條件，得到的統(tǒng)計結果如表1 所示。

表1 海底地形模型與多波束數(shù)據(jù)差值統(tǒng)計結果Table 1 Difference statistics of seafloor topographic model with multibeam data

以標準差作為質量評估依據(jù)，SVG 模型相較于SGA 模型，其精度提高了13.11%，表明在地形起伏變化劇烈的區(qū)域，重力垂直梯度異常能更好的反映地形起伏變化；使用傳統(tǒng)最小二乘賦權法構建的模型MGM_LS，其精度低于SVG 模型，而通過滑動窗口賦權構建的模型MGM 相較于僅使用單一重力數(shù)據(jù)構建的SGA 和SVG 模型，精度分別提高了14.92%和2.08%左右，表明本文采用的滑動窗口賦權融合了兩種重力數(shù)據(jù)在不同地形起伏變化下的反演優(yōu)勢，在復雜海域的水深反演表現(xiàn)更為突出。

4.3 MGM 精度評價

船測水深數(shù)據(jù)在預處理后，可作為研究海域內的真值對反演模型進行外部檢核。將MGM 內插至船測水深處與船測水深值做差，并定義模型與檢核點海深差值與檢核點海深值之比作為相對誤差，得到的差值分布如圖8 所示。

結合圖8a 和圖8b，MGM 與單波束測深數(shù)據(jù)在絕大部分海域差值較小，與多波束測深數(shù)據(jù)差值較大的區(qū)域分布在地形起伏變化劇烈的海山鏈、海隆與海溝交界處；根據(jù)直方圖統(tǒng)計，差值大于1 000 m 的僅占4.48%，可能是由于多波束實測數(shù)據(jù)集中分布在地形變化幅度大的海域，存在粗差未剔除；由圖8c 可知，在兩側相對平坦海域，MGM 相對誤差保持在3%以下，精度較高；在中心海山鏈區(qū)域，相對誤差在10%～50%范圍內占比約為7%，是模型的主要誤差源。

圖8 MGM 差值分布圖Fig.8 Difference distribution diagram of multiple gravity model

為進一步驗證MGM 的可靠性，引入DTU18 和V23.1 模型進行比較分析，得到的統(tǒng)計結果如表2所示。

表2 海底地形模型在檢核點處差值統(tǒng)計結果Table 2 Difference statistics of seafloor topographic model at check point

由表2 知，在研究海域，本文聯(lián)合重力異常和重力垂直梯度異常構建的模型MGM，其誤差極值更小，相較于DTU18 和V23.1 模型，其精度提高了約36%。其次，MGM 與船測數(shù)據(jù)間相關系數(shù)達0.998 0，表明MGM 與實測水深數(shù)據(jù)更為貼近。由于3 種模型中均存在與船測水深數(shù)據(jù)差值較大的點，但僅占檢核點總數(shù)的極小部分，為直觀看出其分布情況，取相對誤差大于6%的點位作為精度較低的異常點。分別統(tǒng)計模型與檢核點差值在400 m 內的頻率分布直方圖（圖9）以及精度較低的異常點位分布圖（圖10）。

圖9 檢核點處差值結果統(tǒng)計直方圖Fig.9 Statistical histogram of the difference results at the check point

圖10 異常點空間位置分布圖（黑色圓點代表異常點位）Fig.10 Spatial distribution map of outliers (black dots represent outliers)

根據(jù)直方圖可看出，MGM 與船測水深數(shù)據(jù)一致性良好，兩者差值在300 m 以內的占比達99.42%，優(yōu)于DTU18 和V23.1 模型的98.12%和98.66%，相較于以上兩種模型，MGM 的異常點數(shù)量分別減少了21.18%和17.32%，表明在地形起伏變化劇烈海域，MGM 仍與船測水深數(shù)據(jù)間有較好的擬合效果。基于前文分析，3 種模型在異常點位的空間位置分布上具有一致性，主要集中在中心海山鏈附近，該部分海域地形起伏變化劇烈，對反演地形容易造成較大的影響。

地形剖面圖可以直觀反映出剖面線上的地形起伏情況。為進一步探究皇帝海山鏈區(qū)域對反演精度的影響，提取其地形剖面進行比較分析[29]，結果如圖11所示。其中藍色曲線代表MGM 繪制的地形剖面圖，紅色曲線代表船測水深數(shù)據(jù)，可作為該地形剖面的實際地形。

圖11 皇帝海山鏈地形剖面圖Fig.11 The topographic profile of the Emperor Seamount Chain

由圖11 知，皇帝海山鏈區(qū)域的整體地形走勢清晰明顯，該剖面的海深范圍為984～6 475 m，平均水深約為4 428 m。MGM 與實測數(shù)據(jù)的差異主要分布在42°～44°N 海域附近，該區(qū)域地形起伏變化劇烈，海山海溝連續(xù)交錯，對反演結果造成了一定程度的影響。整體上，MGM 與實測水深數(shù)據(jù)在描述海山地形走勢上表現(xiàn)出了良好的一致性，兩者相關性達0.996 3。

綜上所述，文中利用滑動窗口賦權聯(lián)合重力異常與重力垂直梯度異常反演海底地形的方法具有較強的可行性，尤其在地形起伏變化劇烈的海域，可有效提高海底地形反演精度。

5 結論

本文利用重力異常和重力垂直梯度異常數(shù)據(jù)，采用滑動窗口賦權和穩(wěn)健回歸分析反演了太平洋皇帝山海域的海底地形模型，將反演結果與船測水深數(shù)據(jù)進行比較分析，得出以下結論：

（1）在海底地形起伏變化劇烈的區(qū)域，依據(jù)重力垂直梯度異常反演的海底地形模型更能準確地反映地形起伏變化。

（2）采用滑動窗口賦權聯(lián)合重力異常和重力垂直梯度異常數(shù)據(jù)，相比于單一數(shù)據(jù)源能有效提高精度，但在地形起伏變化劇烈區(qū)域的細節(jié)表現(xiàn)仍有待提高。

（3）在與船測檢核點比較中發(fā)現(xiàn)，MGM 中精度較低的異常點數(shù)量相較于DTU18 和V23.1 模型有所減少，反映出在地形起伏變化劇烈的區(qū)域，其檢核精度略優(yōu)于DTU18 和V23.1 模型。

采用回歸分析法反演海底地形，對船測數(shù)據(jù)的分布和數(shù)量要求較高，后續(xù)可以針對不同船測數(shù)據(jù)的分布情況以及重力異常和重力垂直梯度異常數(shù)據(jù)在不同地形起伏條件下的適用性開展進一步研究。