范 雕，李姍姍，孟書宇，張金輝，單建晨，王傲明

（1.信息工程大學(xué)，鄭州 450001；2.32022 部隊，武漢 430074）

海洋測繪作為認(rèn)識、開發(fā)、利用海洋的基礎(chǔ)手段和系統(tǒng)工程，測繪成果可為維護(hù)國家海洋權(quán)益、發(fā)展海洋經(jīng)濟、保護(hù)海洋環(huán)境等國家重大發(fā)展戰(zhàn)略順利實施提供基礎(chǔ)數(shù)據(jù)支撐。海底地形測量作為海洋測繪重要組成部分，數(shù)據(jù)及產(chǎn)品廣泛應(yīng)用于國民經(jīng)濟建設(shè)、國防建設(shè)以及科學(xué)研究等方方面面。目前海底地形數(shù)據(jù)獲取方式仍主要通過船舶為載體搭載單波束或者多波束等聲學(xué)儀器進(jìn)行測量，面對占地球表面積大約71%的浩瀚海洋，聲學(xué)測深數(shù)據(jù)覆蓋范圍十分有限，如最新版本的全球海底地形模型SRTM15+V2.0，其中收集的全球聲學(xué)海深測量數(shù)據(jù)在空間格網(wǎng)大小為15"環(huán)境下采樣率也僅為10.84%[1]。有學(xué)者指出，若望采用船基聲學(xué)海底地形測量獲取覆蓋全球的深海水深數(shù)據(jù)，需要持續(xù)投入巨額資金，耗費大約200 船年（ship-year）時間方可完成[2]。因此面對海洋觀測系統(tǒng)化、全球化發(fā)展趨勢，探索低成本、高效率的海底地形數(shù)據(jù)獲取手段顯得尤為重要。

地殼均衡補償理論表明海面重力數(shù)據(jù)與海底地形呈現(xiàn)較強相關(guān)關(guān)系，這種相關(guān)關(guān)系使得利用海面重力信息恢復(fù)海底地形成為可能[3]。特別是隨著測高衛(wèi)星、重力衛(wèi)星等空間大地測量技術(shù)不斷發(fā)展，依托星基平臺高效快速地恢復(fù)了高精度全球海洋重力場信息[4,5]，從而國內(nèi)外學(xué)者對采用海面重力數(shù)據(jù)反演海底地形開展了大量有益研究[6-9]，目前常用的依據(jù)海面重力數(shù)據(jù)反演海底地形算法主要有重力地質(zhì)法（GGM 方法）[10,11]，回歸分析技術(shù)[12-14]，導(dǎo)納函數(shù)方法[15-17]，最小二乘配置技術(shù)[18,19]，模擬退火方法[20]等。其中利用海底地形和海面重力信息線性關(guān)系反演海底地形（稱為回歸分析技術(shù)）的關(guān)鍵核心是確定地形-重力比例因子。關(guān)于比例因子確定問題，國內(nèi)外學(xué)者做了大量的研究探索工作，文獻(xiàn)[12][21][22]采用最小二乘的方法獲取比例因子，文獻(xiàn)[23][24]直接將船測點殘余海底地形和殘余重力數(shù)據(jù)相除得到比例因子，然而若輸入數(shù)據(jù)存在粗差，采用此兩類比例因子獲取方法會嚴(yán)重影響海底地形反演質(zhì)量。文獻(xiàn)[25]文獻(xiàn)為了解決該問題，采用反演波段海底地形和反演波段重力數(shù)據(jù)的中位數(shù)之比獲得比例因子，并強制地形-重力過原點，該方法具有一定的抗差性，但是忽略了地形-重力關(guān)系中常數(shù)項對海底地形模型構(gòu)建結(jié)果的影響。

基于以上分析，兼顧輸入數(shù)據(jù)中粗差對海底地形反演結(jié)果干擾和地形-重力關(guān)系中常數(shù)項對海底地形模型構(gòu)建影響，本文提出采用抗差估計理論確定海底地形-重力比例因子，進(jìn)而構(gòu)建合理海底地形模型的新思路。然后以日本海部分海域作為算法驗證區(qū)，分別利用試驗海區(qū)衛(wèi)星測高重力異常、衛(wèi)星測高重力異常垂直梯度數(shù)據(jù)，結(jié)合稀疏船測水深數(shù)據(jù)，建立了相應(yīng)的海底地形模型，并將本文海底地形模型反演結(jié)果與采用最小二乘方法、中位數(shù)方法構(gòu)建的海底地形模型、DTU10 模型和ETOPO1 模型進(jìn)行了模型橫向比對和精度評估。

1 原理和方法

1.1 基本理論

Parker 公式推導(dǎo)了以海底地形為輸入，海面重力異常為輸出的頻率域表達(dá)式：

式中：G為地球引力常數(shù)；Δρ為海水和地殼的密度差異；d為平均海深；F（Δg(x,y)）和F（h(x,y)）分別表示重力異常和海底地形的二維傅里葉變換；k為波數(shù)kx和ky分別表示x和y方向的波數(shù)（kx=2π/λx，ky=2π/λy）；λx、λy分別表示x和y方向的波長。當(dāng)僅考慮式(1)中線性項時，海底地形和海面重力異常的頻域線性關(guān)系可近似表示為

式中，Z(kx,ky)表征將海底地形轉(zhuǎn)換為重力異常的能力，通常稱為重力導(dǎo)納函數(shù)。重力導(dǎo)納函數(shù)Z(kx,ky)具有各向同性、空間不變性特點。研究表明，大尺度的海底地形特征可以用地殼均衡補償理論來解釋，那么當(dāng)顧及地殼均衡補償影響時，式(2)可表示為

式中，Zf(kx,ky)稱為撓曲導(dǎo)納函數(shù)；φ(kx,ky)為均衡響應(yīng)函數(shù)：

式中，Tc為地殼厚度；Φe(kx,ky)為撓曲均衡響應(yīng)函數(shù)（Flexural Response Function）。分析(3)式可知，均衡響應(yīng)函數(shù)將由地形產(chǎn)生的重力異常修正為地形和均衡補償物質(zhì)對重力異常的貢獻(xiàn)之和。撓曲均衡響應(yīng)函數(shù)定義為：

其中，g為研究區(qū)域重力值；D為撓曲剛度，計算式如下：

式中，E為楊氏模量（Young’s Modulus）；Te為板塊的有效彈性厚度；υ為泊松比。據(jù)式(3)可得利用重力異常反演海底地形算法模型為

式中，轉(zhuǎn)換函數(shù)Qf(kx,ky)為導(dǎo)納函數(shù)Zf(kx,ky)的倒數(shù)。為了抑制均衡補償對海面重力異常的影響和克服重力異常向下延拓過程（轉(zhuǎn)換函數(shù)Qf(kx,ky)包含向下延拓因子ekd）噪聲放大的缺點以獲得穩(wěn)定的海底地形反演結(jié)果，引入帶通濾波器W(kx,ky)限制轉(zhuǎn)換函數(shù)Qf(kx,ky)，則

式中，h0為對應(yīng)波段的海底地形。那么在未考慮地殼均衡補償環(huán)境下，同時根據(jù)頻率域求導(dǎo)法則，可得重力異常/重力異常垂直梯度（Δgz=-(?Δg/ ?z)）為系統(tǒng)輸入，海底地形為系統(tǒng)輸出的頻率域表達(dá)式

由式(9)可知，重力異常/重力異常垂直梯度經(jīng)濾波等處理后與海底地形呈線性比例關(guān)系，比例因子為S(x,y)=(1/2π GΔρ)。從而可以認(rèn)為海底地形反演工作通常在中長波波段范圍開展，而短波和長波部分海底地形信息可通過船測海深等方式進(jìn)行補充，進(jìn)而獲得全波段的海底地形

式中，hl(x,y)為非反演波段海底地形；Δg0(x,y)為經(jīng)濾波且向下延拓后的重力異常；Δgz0(x,y)為經(jīng)濾波且向下延拓等處理后的重力異常垂直梯度相關(guān)量，其他參數(shù)意義同上。經(jīng)分析可知，依據(jù)式(10)反演海底地形，關(guān)鍵在于確定比例因子。然而實際上由于海底地形地貌構(gòu)造復(fù)雜和沉積層密度異常等因素影響，比例因子并非常數(shù)而是一個變化量[25]，不能根據(jù)理論密度參數(shù)直接計算，而是通過船測數(shù)據(jù)約束，分別計算不同區(qū)域的線性參數(shù)[22]。

圖1 海底地形反演流程圖Fig.1 Flow chart of seafloor topography inversion

基于以上分析，繪制海底地形的反演流程如圖1所示。需要說明的是，傳統(tǒng)方法首先獲取控制點處離散比例因子，然后將比例因子格網(wǎng)化處理，進(jìn)而構(gòu)建海底地形模型[25]。文獻(xiàn)[12]指出，依據(jù)移動窗口法獲取比例因子構(gòu)建海底地形質(zhì)量高于傳統(tǒng)方法。據(jù)此，本文采用移動窗口法求取比例因子參數(shù)。

1.2 比例因子獲取方法

1.2.1 中位數(shù)方法（Median 方法）

假設(shè)反演海區(qū)有限頻段海底地形和重力相關(guān)數(shù)據(jù)分別為b0(x,y) 和g0(x,y)，則二者比例因子可表示為：

式中，med(b0(x,y))和med(g0(x,y))分別表示地形和重力相關(guān)數(shù)據(jù)的中位數(shù)。

1.2.2 最小二乘方法（Ls 方法）

為獲得海底地形-重力的比例因子，以重力數(shù)據(jù)為自變量，海底地形為因變量進(jìn)行線性擬合獲取比例因子和常數(shù)項參數(shù)結(jié)果。依據(jù)最小二乘方法解算參數(shù)結(jié)果為

1.2.3 抗差估計方法（Ro 方法）

抗差估計（RobustEstimation）也稱為穩(wěn)健估計或者魯棒估計，源于統(tǒng)計學(xué)中的抗差性（Robustness）概念。當(dāng)觀測值中存在粗差的情況下，采用抗差估計可以得出正常模式下的最佳估值?？共罟烙嫷脑瓌t是充分利用有效信息，限制利用有用信息，排除有害信息，在所假定的模型下，可以獲得可靠、有效、具有實際意義的參數(shù)估值。

依據(jù)抗差估計理論，參數(shù)解算迭代策略為

式中，k為迭代次數(shù)；為擬合參數(shù)結(jié)果；A為重力數(shù)據(jù)構(gòu)成的系數(shù)矩陣；L為海底地形輸入向量；為等價權(quán)矩陣。中的元素為

式中，稱為等價權(quán)；ω(vi)稱為權(quán)因子；ψ(vi)是非線性函數(shù)；vik表示第k次迭代的殘差分量。實際計算等價權(quán)的操作中，首先需將殘差vi標(biāo)準(zhǔn)化

式中，miv是vi的中誤差；Ai表示系數(shù)陣A的第i行向量；σ0是單位權(quán)中誤差，采用標(biāo)準(zhǔn)化殘差絕對值的中位數(shù)計算

當(dāng)參數(shù)x 的迭代解滿足式(17)，迭代停止，獲得最終解。

式中，ε為設(shè)置的迭代停止條件。本文選擇的等價權(quán)函數(shù)為Huber 函數(shù)

由式(18)可知，當(dāng)≤2σ0，Huber 函數(shù)等價于最小二乘方法，當(dāng)＞2σ0，Huber 函數(shù)利用L1范數(shù)對觀測值進(jìn)行降權(quán)，這樣可以有效抑制粗差的不利影響，同時保持良好的正態(tài)有效性。

值得說明的是，Ro 方法求解比例因子本質(zhì)核心仍可認(rèn)為是最小二乘估計。不同的是，Ls 方法獲取比例因子僅僅是有船測海深數(shù)據(jù)約束的最小二乘方法；Median 方法使用排序信息獲取比例因子；Ro 方法獲取比例因子是具有外部Huber 約束地形信息的正則化最小二乘法。從而算法本質(zhì)上是希望對最小二乘施加約束，進(jìn)而獲得更有效、合理的比例因子。

2 數(shù)值分析試驗

2.1 數(shù)據(jù)準(zhǔn)備和前期處理

選取日本海4 °×4 °（131 °E～135 °E,36 °N～40 °N）海域作為試驗海區(qū)開展數(shù)值分析試驗。該海域西北方向地形平坦，東北和東南方向存在海山等地形復(fù)雜地貌，從而認(rèn)為該海域具備作為算法可靠性驗證的典型區(qū)域條件。試驗過程以船測點作為約束條件補充海底地形短波部分和長波部分地形信息。數(shù)據(jù)來源如下：①船測海深數(shù)據(jù)：來源于 NGDC（The National Geophysical Data Center）發(fā)布的研究海區(qū)實測水深數(shù)據(jù)。研究海域原始船測水深數(shù)據(jù)共39018 個，通過3σ準(zhǔn)則對實測水深數(shù)據(jù)進(jìn)行粗差剔除，最終得到37750個實測水深數(shù)據(jù)。選取其中約五分之四的數(shù)據(jù)（30877個）作為控制點，以期對構(gòu)建的海底地形模型進(jìn)行控制；剩余6873 個實測水深數(shù)據(jù)作為外部檢核點供海底地形模型精度評估使用?？刂泣c和檢核點的分布如圖2(a)所示，其中白色六邊形為控制點，黑色三角形為檢核點，圖中背景為ETOPO1 海深模型。②衛(wèi)星測高重力數(shù)據(jù)包括重力異常和重力異常垂直梯度數(shù)據(jù)均來自于SIO，UCSD（Scripps Institution of Oceanography,University of California,San Diego），V24.1 版本，分辨率為1’，研究海域重力異常和重力異常垂直梯度數(shù)據(jù)如圖2(b)(c)所示。

圖2 研究海域數(shù)據(jù)Fig.2 Study sea area data

2.2 海底地形模型構(gòu)建

根據(jù)利用衛(wèi)星測高重力數(shù)據(jù)反演海底地形理論可知，海底地形模型構(gòu)建過程中，由于短波部分重力數(shù)據(jù)所含海底地形信噪比較低以及重力數(shù)據(jù)向下延拓存在噪聲放大現(xiàn)象，從而采用低通濾波器對重力數(shù)據(jù)進(jìn)行濾波處理以消除短波部分反演結(jié)果不穩(wěn)定影響。為盡可能使用重力信息，同時考慮低通濾波器和轉(zhuǎn)換函數(shù)特點以及反演分辨率等因素，結(jié)合文獻(xiàn)[21]和文獻(xiàn)[25]，最終設(shè)定短波截止波長為15 km。另外，為有效抑制地殼均衡補償對海面重力信息的影響，充分利用衛(wèi)星測高重力信息中的地形分量，依據(jù)文獻(xiàn)[26]認(rèn)為研究海域有效彈性厚度約為5 km，然后根據(jù)有效彈性厚度與波長的關(guān)系，長波截止波長最終設(shè)定為160 km。基于以上分析，日本海試驗海域利用衛(wèi)星測高重力數(shù)據(jù)反演海底地形的反演波段范圍為15 km～160 km。

分別以重力異常和重力異常垂直梯度作為系統(tǒng)輸入，采用Median 方法、Ls 方法和Ro 方法計算試驗海區(qū)線性分析一次項和常數(shù)項參數(shù)結(jié)果分別如圖3～圖4所示。

圖3 試驗海區(qū)topo-vgg 線性分析結(jié)果（無粗差）Fig.3 Topo-vgg linear analysis results in experimental seaarea (no gross error)

圖4 試驗海區(qū)topo-dg 線性分析結(jié)果（無粗差）Fig.4 Topo-dg linear analysis results in experimental sea area(no gross error)

最終，以重力異常和重力異常垂直梯度為輸入源，依據(jù)各方法獲得的回歸參數(shù)結(jié)果構(gòu)建相應(yīng)的海底地形模型如圖5和圖6所示。

觀察圖5和圖6海深模型構(gòu)建結(jié)果，六種海深模型均展現(xiàn)了研究海域地形地貌特點：海域西北部以海底平原為主要地形呈現(xiàn)；海域東北和東南部地貌形態(tài)豐富，其中海山大面積分布于該區(qū)域，海域東部邊緣存在海底平原地貌。仔細(xì)對比圖5和圖6海深模型可以發(fā)現(xiàn)，使用Median 方法，以重力異常為輸入數(shù)據(jù)構(gòu)建的海深模型細(xì)節(jié)呈現(xiàn)明顯豐富，如圖5(a)西北部可以清晰看見3 處地形隆起，而圖5(b)(c)、圖6(a)～(c)目視僅1 處地形隆起，另外兩處則較為模糊；又如在（39 ° N～40 ° N，133 ° E～134 ° E）區(qū)域，圖5(a)中地形呈現(xiàn)細(xì)節(jié)明顯好于其余模型。

圖5 重力異常構(gòu)建海深模型Fig.5 Bathymetry model derived by gravity anomaly

圖6 重力異常垂直梯度構(gòu)建海深模型Fig.6 Bathymetry model derived by vertical gravity gradient anomaly

分析沉積層對比例因子和相關(guān)系數(shù)結(jié)果影響，依據(jù)CRUST 1.0 繪制試驗海域沉積層厚度和密度分布圖如圖7所示。比對圖3、圖4和圖7可大致發(fā)現(xiàn)，海底地形-重力比為適應(yīng)沉積層厚度和密度變化而變化：在沉積層厚度較薄區(qū)域，沉積層密度較小，海水和地殼密度差異較大，從而海底地形-重力呈現(xiàn)出較強相關(guān)性，如研究海域東北部區(qū)域；沉積層厚度較厚區(qū)域，沉積層密度較大，海底地形起伏較小，覆蓋在地形表面的沉積層也會對海面重力異常產(chǎn)生影響，該區(qū)域地形-重力相關(guān)性較弱，如西南部和西北部研究海域。另外，對于海底地形和重力數(shù)據(jù)比例因子較大的原因，筆者認(rèn)為可能是，因地殼均衡補償作用對海面重力影響，導(dǎo)致該海區(qū)重力值（重力異常/重力異常垂直梯度）較小，從而海底地形和重力異常比例因子較大。雖然實驗開始階段輸入數(shù)據(jù)經(jīng)過濾波處理，然而因地形撓曲波長的選擇為試驗海區(qū)大范圍內(nèi)的平均值，當(dāng)部分海域地形波長可能小于撓曲波長時，該區(qū)域可能發(fā)生局部均衡補償現(xiàn)象，從而導(dǎo)致海面重力受到均衡補償影響。此時可通過縮小濾波窗口，減小撓曲波長的方法抑制均衡補償影響，然而縮小濾波窗口將引發(fā)其他問題。

圖7 試驗海域沉積層Fig.7 Sediment in study area

2.3 精度評價分析

以檢核點作為外部檢核條件，引入國際廣泛使用的ETOPO1 海深模型和DTU10 海深模型進(jìn)行比較，評估各海深模型精度和效能。采用雙線性插值方法將模型海深值內(nèi)插到外部檢核點處，并與檢核點處實測水深值作差以評估海深模型精度。

統(tǒng)計結(jié)果表明，以重力異常為數(shù)據(jù)源，采用Median 方法建立的BGMed 海深模型檢核差值最大值和最小值絕對值均超過 1000 m，檢核均方差為243.37 m，檢核精度遠(yuǎn)低于基于Ls 方法和Ro 方法構(gòu)建的BGLs 模型和BGRo 模型。其中采用Ro 方法反演得到的BGRo 海深模型的檢核均方差為90.03 m，Ls 方法反演海深模型的檢核均方差為94.04 m，可以認(rèn)為利用Ro 方法與Ls 方法在研究海域構(gòu)建的海底地形模型整體效果相當(dāng)，因此，從整體上看，Ro 方法具有良好的適用性。ETOPO1 海深模型和DTU10 海深模型的檢核均方差為99.94 m 和125.42 m，檢核精度高于BGMed 模型，低于BGLs 模型和BGRo 模型。采用重力異常垂直梯度為數(shù)據(jù)源建立的BVGMed 模型、BVGLs 模型和BVGRo 模型中，BVGRo 海深模型檢核精度最高，為89.51 m，BVGMed 海深模型檢核精度最低，為264.77 m。BVGRo 模型相較于ETOPO1海深模型和DTU10 海深模型精度分別提高了10.44%和28.63%左右。對比不同衛(wèi)星測高重力數(shù)據(jù)構(gòu)建的海深模型，BVGLs 海深模型精度優(yōu)于BGLs 海深模型，BGRo 海深模型檢核精度低于BVGRo 海深模型。由于重力異常和重力異常垂直梯度隨距離的衰減特性不同，為探究二者面對淺海域和深海域表現(xiàn)狀態(tài)，繪制BGRo 檢核差值絕對值和BVGRo 檢核差值絕對值之差與海深關(guān)系如圖8所示?？疾靾D8可以發(fā)現(xiàn)，本文試驗區(qū)域水深小于1500 m 左右深度，BVGRo 檢核差值絕對值比較明顯小于BGRo 檢核差值絕對值，即該水深范圍內(nèi)重力異常垂直梯度反演效果較好；水深大于1500 m 左右范圍，BGRo 檢核差值絕對值比較明顯小于BVGRo 檢核差值絕對值，即該水深范圍內(nèi)重力異常反演效果較佳。筆者認(rèn)為從重力異常和重力異常垂直梯度數(shù)據(jù)的特性角度考慮：重力異常反映的是相對低頻的信息，而重力異常垂直梯度反映的是相對高頻的信息，因此海域較深，重力異常數(shù)據(jù)的反演結(jié)果相對較好，而淺部則是重力異常垂直梯度較好。

就外部檢核結(jié)果統(tǒng)計指標(biāo)整體而言，Ro 方法和Ls 方法在試驗海域構(gòu)建模型精度相當(dāng)，優(yōu)于Median方法；而Median 方法對于試驗海域地形呈現(xiàn)豐富度又好于Ro 方法和Ls 方法?？梢?，三種方法各自具有不同的海底地形構(gòu)建優(yōu)勢。

圖8 BGRo 檢核差絕對值和BVGRo 檢核差絕對值之差與海深關(guān)系Fig.8 Relationship between the difference between checking difference absolute value of BGRo and BVGRo and sea depth

分別統(tǒng)計以重力異常和重力異常垂直梯度為數(shù)據(jù)源，采用Median 方法、Ls 方法和Ro 方法構(gòu)建的海深模型不同檢核差值范圍內(nèi)檢核點數(shù)量，結(jié)果如圖9所示。圖9顯示Median 方法構(gòu)建的海底地形模型檢核差值在不同差值區(qū)間數(shù)量穩(wěn)定；Ro 方法構(gòu)建的海深模型檢核差值小于50 m 左右范圍內(nèi)檢核點數(shù)量明顯多于Ls 方法和Median 方法，檢核差值大于50 m 范圍內(nèi)檢核點數(shù)量明顯少于Ls 方法和Median 方法，進(jìn)一步說明了Ro 方法良好的估計性質(zhì)和適用性。進(jìn)一步驗證Ro 方法的優(yōu)越性，首先以試驗海域反演波段海底地形和重力異常/重力異常垂直梯度為研究對象，分別選取其中3025 個數(shù)據(jù)點進(jìn)行分析研究，包括反演波段海底地形、反演波段重力異常和重力異常垂直梯度數(shù)據(jù)。試驗過程中認(rèn)為原始反演波段輸入數(shù)據(jù)中不含粗差，原始反演波段海底地形、反演波段重力異常和重力異常垂直梯度輸入數(shù)據(jù)（分別簡稱為topo、dg 和vgg）統(tǒng)計結(jié)果如表1。實驗過程中認(rèn)為重力數(shù)據(jù)質(zhì)量良好無粗差存在[25]，分別在與dg 和vgg相擬合的topo 中隨機加入粗差（在原反演波段地形數(shù)據(jù)基礎(chǔ)上增加或減小5～10 倍中誤差作為粗差），最終與dg 相擬合的topo 數(shù)據(jù)中加入粗差點數(shù)為24 個，約占總點數(shù)的7.93‰，得到含粗差的地形數(shù)據(jù)topoErrDg；與vgg 相擬合的topo 數(shù)據(jù)加入粗差的點數(shù)為12 個，約占總點數(shù)的 3.97‰，獲得含粗差的地形數(shù)據(jù)topoErrVgg。添加粗差的地形數(shù)據(jù)統(tǒng)計結(jié)果見表1。相較于添加粗差前的topo 數(shù)據(jù)，添加粗差后的地形數(shù)據(jù)最值變化明顯，數(shù)據(jù)整體波動劇烈。

圖9 模型比較Fig.9 Bathymetry model comparison

表1 輸入數(shù)據(jù)統(tǒng)計結(jié)果Tab.1 Statistical results of input data

分別采用Median 方法、Ls 方法和Ro 方法確定地形-重力比例因子，其中反演波段海底地形和重力異常/重力異常垂直梯度線性回歸分別簡稱為topo-dg 和topo-vgg。未加粗差的情況下，采用三種線性分析技術(shù)獲得海底地形和重力異常/重力異常垂直梯度的擬合結(jié)果如圖10(a)(c)所示，圖10(a)(c)顯示海底地形和重力表現(xiàn)良好的線性性質(zhì)。在地形數(shù)據(jù)隨機加入粗差后與相應(yīng)的重力數(shù)據(jù)（dg 和vgg）采用相同的方法獲得的線性分析結(jié)果如圖10(b)(d)所示，圖10(b)(d)中藍(lán)色圓點表示加入粗差后的地形數(shù)據(jù)結(jié)果。圖10(b)(d)顯示，加入粗差對抗差估計結(jié)果影響較小。

圖10 海底地形和重力數(shù)據(jù)Fig.10 Seafloor topography and gravity data

三種方法擬合參數(shù)統(tǒng)計結(jié)果表明，反演波段海底地形數(shù)據(jù)加入粗差后，海底地形和重力數(shù)據(jù)相關(guān)性明顯降低，如未加粗差的海底地形和重力異常垂直梯度的相關(guān)系數(shù)為0.76，加入粗差的海底地形和重力異常垂直梯度相關(guān)系數(shù)為0.65。因此海底地形和重力回歸分析過程中應(yīng)充分利用有效信息，限制利用有用信息，排除有害信息，進(jìn)而盡量避免粗差對結(jié)果的影響。經(jīng)計算，不含粗差情況下，Median 法、Ls 法和Ro 法獲得的比例因子分別為28.54、14.95 和16.29；加入粗差后三種方法獲得的比例因子分別為28.59、13.78 和16.32，與不含粗差環(huán)境相比，分別在原值基礎(chǔ)上變動1.75‰、7.83%和1.84‰，由此可見粗差對最小二乘影響較大，對抗差估計結(jié)果影響較小。同理可得加入粗差后，三種線性分析方法解算的海底地形和重力異常垂直梯度的比例因子較不含粗差情況分別變動0.00%、2.42%和0.00%，體現(xiàn)了抗差估計在粗差環(huán)境中優(yōu)良的線性分析性質(zhì)。

基于以上分析，在每個移動窗口內(nèi)的初始反演波段海底地形輸入數(shù)據(jù)中隨機加入粗差（粗差點占總數(shù)的3%～6%），然后分別依據(jù)Median 方法、Ls 方法和Ro 方法構(gòu)建輸入海底地形數(shù)據(jù)含粗差的海深模型，構(gòu)建的海深模型檢核統(tǒng)計結(jié)果如表2。

表2 海深模型（含粗差）差值檢核統(tǒng)計結(jié)果（單位：米）Tab.2 Statistical checking results of difference bathymetry model (including gross errors) (Unit:m)

為描述方便，以重力異常為輸入源，采用Median方法、Ls 方法和Ro 方法建立的海深模型分別稱為BGMedErr 模型、BGLsErr 模型和BGRoErr 模型；以重力異常垂直梯度為輸入數(shù)據(jù)，采用Median 方法、Ls 方法和 Ro 方法構(gòu)建的海深模型分別稱為BVGMedErr 模型、BVGLsErr 模型和BVGRoErr 模型。表2檢核統(tǒng)計結(jié)果顯示，在輸入海底地形數(shù)據(jù)中添加粗差，將不同程度地影響海深模型檢核精度。如以重力異常為系統(tǒng)輸入，BGMedErr 模型檢核差值均方差為244.21 m，BGMed 模型與檢核點差值均方差為243.37 m，反映了Median 方法具有一定的抵抗粗差能力。但是采用Median 方法獲得的海深模型，由于反演過程強制擬合過原點，忽略了常數(shù)項對結(jié)果的影響，最終的海深模型反演精度低于 BGLsErr 模型和BGRoErr 模型。BGLsErr 模型相較于BGLs 模型檢核均方差變化了16.88 m，而BGRoErr 模型相較于BGRo模型檢核均方差僅變化了0.02 m，證明了采用Ro 方法反演海深模型能夠充分利用有效信息，排除有害信息，獲得穩(wěn)定、有效的解算結(jié)果。以重力異常垂直梯度構(gòu)建的海深模型檢核統(tǒng)計結(jié)果也表明，Median 方法獲得海深模型精度最低，Ro 方法獲取的海深模型檢核精度最高。同時BVGLsErr 模型相較于BVGLs 模型檢核均方差變化了11.73 m，而BVGRoErr 模型相較于BVGRo 模型檢核均方差僅變化了0.25 m。

各個海深模型加入粗差前后不同檢核差值區(qū)間檢核點變化情況如圖11所示?？砂l(fā)現(xiàn)加入粗差后對Ls 方法構(gòu)建的海深模型影響相對較大：不同檢核差值區(qū)間檢核點數(shù)量變化明顯，檢核差值小的檢核點數(shù)量減小，檢核差值大的檢核點數(shù)量增多；而加入粗差后Ro 方法構(gòu)建的海深模型不同檢核差值區(qū)間檢核點數(shù)量并未發(fā)生較大變化，進(jìn)一步證明了Ro 方法對待含粗差輸入數(shù)據(jù)具備較高穩(wěn)定性。

圖11 加入粗差后模型比較Fig.11 Bathymetry model comparison after adding gross errors

3 結(jié) 論

針對當(dāng)前利用回歸分析技術(shù)恢復(fù)海底地形算法中有限波段重力數(shù)據(jù)與海底地形之間比例因子獲取忽略常數(shù)項或者抵抗粗差能力弱缺陷，提出了可采用抗差估計理論，實現(xiàn)兼顧輸入數(shù)據(jù)粗差影響和常數(shù)項的比例因子參數(shù)求解新思路。選取日本海部分海域作為算法驗證區(qū)，以衛(wèi)星測高重力異常和重力異常垂直梯度為重力輸入，分別使用最小二乘法、中位數(shù)法和抗差估計法恢復(fù)了試驗海域海底地形；同時討論了初始海底地形存在粗差環(huán)境下，三種比例因子獲取方法抵抗粗差影響能力強弱，引入ETOPO1 和DTU10 海深模型對本文模型進(jìn)行了橫向比對分析。結(jié)果表明，由于重力異常垂直梯度在較高頻段所含“地形”信息相比重力異常更加豐富，從而研究海域以重力異常垂直梯度作為輸入源構(gòu)建的海底地形模型在淺海域整體質(zhì)量優(yōu)于以重力異常作為輸入源構(gòu)建的海底地形模型；同理，深海域重力異常建立海底地形效果較好。中位數(shù)方法、最小二乘方法和抗差估計方法構(gòu)建的試驗海域海底地形模型中，由于中位數(shù)方法強制回歸分析過原點，忽略常數(shù)項的影響，檢核精度最低；抗差估計構(gòu)建的海底地形模型檢核精度最高。最小二乘方法和抗差估計方法構(gòu)建海底地形模型地貌呈現(xiàn)整體效果相當(dāng)，而中位數(shù)方法構(gòu)建的海底地形模型地貌細(xì)節(jié)呈現(xiàn)更為豐富。輸入海底地形數(shù)據(jù)含粗差的情況下，抗差估計方法構(gòu)建的海底地形模型檢核精度變化最小，反映了利用抗差估計構(gòu)建海底地形模型具有可靠性高、實用性強特點。

綜上，依托回歸分析方法恢復(fù)試驗海域海底地形時，建議可根據(jù)實際需求，重點考慮使用重力異常垂直梯度，采用抗差估計方法構(gòu)建海底地形模型。