基于最小二乘配置的陌生海域海底地形反演方法

范 雕,李姍姍,馮進凱,黃 炎,范昊鵬,張金輝,李新星

1. 信息工程大學,河南 鄭州 450001; 2. 軍事科學院國防科技創新研究院,北京 100071

海洋是人類生存和發展的重要資源寶庫,其面積占地球表面積71%左右[1]。文獻[2]研究表明目前依據船載聲學測量手段僅有15%左右的海域面積完成了空間分辨率約800 m的海底地形繪制,人類對于海洋的認識仍處于較低水平和層次[3]。文獻[4]認為依托船舶為載體實施海底地形測量,若以空間網格大小為1′作為考慮對象,全球大約還有90%的海域屬于未被勘探狀態,而且船舶測線間存在巨大的空間間隔,測量結果分布很不均勻,人類對海底地形的描繪詳細程度遠不及對月球等其他行星的了解深入[5]。

然而,海底地形地貌數據作為開發和利用海洋的基礎性地理空間信息數據,廣泛影響其他科學研究和應用實踐。一方面,海底地形變化直接影響海洋環流和混合,從而間接影響地球氣候、生物多樣性和糧食產量。另一方面,海底地形和海底構造的精細化表達有利于研究板塊構造規律,有利于揭示海盆和海底地質運動歷史,有利于預判發生地震、火山爆發和海嘯等區域位置。因此,過去100多年來國內外學者持續不斷努力研究繪制全球范圍海底地形圖[6-7],囿于技術裝備發展水平制約,前期研究進展較為緩慢。直到19世紀90年代中期,由于GEOSAT和ERS-1等執行大地測量任務的測高衛星數據不斷積累和應用,全球海洋區域重力場精度和分辨率不斷提升,依托海洋重力信息反演海底地形研究不斷推進[8-10]。至此,運用空天探測技術獲取全球范圍的海底地形信息方才取得重要突破。文獻[11]利用衛星測高重力異常數據,依據回歸方法(又稱S&S方法)建立了范圍涵蓋南北緯72°之間,分辨率約15 km的海深格網模型。隨后數十年時間,由于模型構建數據量不斷豐富,數據精度不斷提升,學者們在更高精度、更高分辨率及更廣覆蓋范圍的海底地形建模方面進行了積極探索,比如:①豐富的船測海深數據可用于校準和修正預測海深結果[12],以使構建的海深模型網格密集程度更高;②波形重跟蹤技術的廣泛使用,在算法層面不斷提高測高結果精度,從而不斷提升海洋重力場質量[13]。目前利用重力數據反演海底地形算法大致可分為統計法和解析法兩大類:統計法包括回歸分析方法[14]、模擬退火方法[15]、最小二乘配置[16]及神經網絡[17]等;解析方法包括重力地質方法[18-19]、導納函數方法[20]及顧及地形高次項影響的解析反演方法[21-22]等。然而,重力地質方法、回歸分析方法、模擬退火方法在試驗海區需要外部船測水深數據支撐,在無船測水深數據支撐的陌生海區該方法將不再適用,進而極大削弱算法應用范圍?；赑arker公式的解析反演方法不需要外部船測水深數據為補償,文獻[23]利用此方法開展了無船測水深數據(利用先驗海深模型代替)的海底地形建模,但是獲取該方法涉及的諸多地球物理參數并非易事,若全球海域均采用參數經驗值代替,毋庸置疑會造成最終反演精度損失,進而導致反演結果可信度降低。

基于此,為擴展無船測水深數據支持環境下陌生海域海底地形恢復方法,本文提出以衛星測高重力數據和先驗海底地形模型為輸入,運用最小二乘配置方法(統計法),解決無船測數據支撐環境下高效構建陌生海區海底地形模型難題的新思路?；静呗允?首先,在船測水深數據和重力數據已知海區統計解算出反映海區地形的相關協方差模型;然后,采用最小二乘配置方法,以海區重力數據為輸入,構建相應的海底地形模型并檢驗算法適用性;最后,將已知海區數據統計結果(不同類型參量的協方差函數)遷移、推廣并應用于與已知海區相似的船測數據陌生海區,實現僅依靠重力數據的海底地形模型構建。同時將本文反演結果與導納函數反演海底地形結果、DTU18、ETOPO1模型進行比對分析,簡要評估模型精度和效能。

1 原理與方法

1.1 最小二乘配置方法

利用解析方法(例如導納函數)反演海底地形時,算法模型涉及的地球參數求解并非易事。文獻[24]認為海底地形和重力信息為兩類隨機變量,從而依據重力異常數據,采用統計方法恢復海底地形模型為

(1)

式中,b為海底地形;g為海面重力異常信息;μb和μg表示海底地形和海面重力異常的數學期望(可用相應的平均值代替);Cbg為海底地形和重力異常的互協方差矩陣;Cgg為海面重力異常的自協方差。式(1)化簡為

(2)

由式(2)可知,使用海面重力異常數據,運用最小二乘配置方法恢復海底地形信息關鍵在于獲取地形-重力異?；f方差函數和重力異常信息間自協方差函數。其中重力信息的自協方差函數求解在物理大地測量重力擬合推估研究較多[26],目前已有成熟的解決手段。然而海底地形和重力異常分屬于兩種不同類型的參數,因此如何精確求解地形-重力異?；f方差函數成了利用式(2)恢復海底地形的主要難題。

1.2 地形-重力協方差估計

將地形和重力異常作為一維隨機變量參考時,地形-重力異?；f方差的無偏估計通過遍歷數據得到[27]

(3)

(4)

式(3)離散化的結果為

(5)

式中,l表示重力異常和地形間距離,通常取0,1,…,N-1。易知,式(5)重力-地形協方差估計是無偏估計

(6)

式中,E(·)表示期望算子。同理,二維重力-地形互協方差函數無偏估計表示為

(7)

借鑒重力異常擬合推估自協方差函數求解策略[26, 28],對于方形格網研究區域(N行和M列)重力-地形的互協方差可采用式(8)統計縱向和橫向互協方差

(8)

式中,q和qj分別表示橫向距離和縱向距離。最終將所有具有距離l的縱向協方差和橫向協方差取平均值作為協方差Cbg(l)的統計數值

(9)

協方差統計數值結果反映了海區相關數據變化的內在規律,通常情況下需構建強相關函數模擬海區數據協方差統計結果,進而依托強相關函數擬合結果構造相應的(互)協方差矩陣,從而據最小二乘配置算法模型恢復海底地形信息。強相關函數如Inverse協方差函數、Gauss協方差函數、Hirvonen協方差函數、似高斯協方差函數、周期指數函數、Tscherning/Rapp協方差函數和Mortiz協方差函數等[29]。本文選擇較為簡單的似高斯協方差函數作為海區協方差統計備選擬合函數,表達式為

Cbg(l)=C0·e-al

(10)

式中,C0是兩點距離為零對應的協方差函數值,表示協方差函數的變化幅度;l為兩點之間距離;a為常數。式(10)通過恒等變換后為

(11)

從而依據最小二乘方法可求解似高斯協方差函數待求參量a。借鑒以上分析結果,同理可求得重力-地形協方差Cgb和重力異常的自協方差Cgg,從而可獲取以海面重力信息為輸入,采用最小二乘配置方法(式(2))構建相應波段的海底地形模型。

2 數值分析試驗

2.1 數據準備與前期處理

選擇西太平洋某1°×1°(155.2°E—156.2°E,22.5°E—23.5°E)海區作為研究試算區域,海區輸入數據包括衛星測高重力異常數據和稀疏船測數據,數據來源和預處理如下。

(1) 海區重力異常來源于SIO(Scripps Institution of Oceanography)于2022年8月3日發布的S&S V32.1重力異常模型,該版本重力異常模型增加了Altika、Cryosat LRM、Cryosat SAR和Sentinel-3A/B等數據。研究海區重力異?？臻g分布如圖1所示,相應的統計結果見表1。重力異常較大值集中在海區中部區域,推測該區域可能存在明顯的海山地形。

表1 研究海域重力異常(S&S V28.1)數據統計

圖1 研究海區重力異常Fig.1 Gravity anomaly in study sea area

(2) 稀疏船測海深數據來源于NGDC(The National Geophysical Data Center)發布的單波束和多波束船測海深數據。 共收集到1061個實測船測點,數據頻數分布直方圖如圖2(a)所示。依圖2(a)可知,研究海區水深測量結果主要集中于5000 m以上。采用一定的處理手段對原始水深測量數據進行質量控制,本文處理策略是依據2019年12月發布的海洋空間格網大小為15″的SRTM15+V2.0作為參考模型[30],采用雙線性內插技術將SRTM15+V2.0模型海深值內插到離散船測點處,并與船測點海深值作差得到海深殘差;借鑒3σ誤差處理方法,認為海深殘差絕對值大于3倍海深殘差標準差的船測點是粗差點,然后作剔除處理。最終利用3σ方法處理原始船測海深數據后,共剔除34個可疑水深點,剩余1027個船測水深點。原始水深、處理后水深及水深殘差統計結果見表2,表中括號內數值代表參與統計的水深點個數。由表2可知,以SRTM15+V2.0作為對比模型,海深殘差標準差為44.59 m,從而粗差閾值為133.77 m。原始船測水深數據剔除粗差后結果的相應直方圖如圖2(b)所示,船測數據空間分布情況如圖3所示(圖中背景為S&S V19.1海深模型)。

表2 海深數據統計結果

圖2 船測數據分布直方圖Fig.2 Distribution histogram of shipborne bathymetry data

圖3 研究海區粗差剔除后船測點分布Fig.3 Distribution of survey points after gross error elimination

利用衛星測高重力數據通常僅可恢復有限波段海底地形信息,剩余頻段地形信息常使用實際船測水深數據進行補償; 另外為評估最終構建的海深模型精度,通常選擇部分船測點作為外部檢核參考[21]。從而本文基于1027個粗差處理后的船測點,隨機均勻選擇847個船測點(約占船測點總數的82.47%)作為控制點以補償非反演頻段地形信息;剩余180個船測點(約占船測點總數的17.53%)作為評估海深模型的外部檢核點。

另外,S&S V19.1海深模型中奇數海深值由實際海深約束生成,偶數海深值為重力反演預測水深值,從而可認為這些奇數海深值是實際船測海深值。研究海區S&S V19.1海深模型中奇數(實際船測數據約束)和偶數(重力預測海深)海深值空間分布情況如圖4(a)所示,其中“·”代表重力預測海深的位置,“?”代表實際水深約束生成的海深值位置。S&S V19.1模型在試驗海區共有616個船測實際值,相應的統計結果見表2,海深平均值和標準差分別為4 861.00 m和816.09 m;船測數據頻數分布直方圖如圖4(b)所示。

圖4 S&S V19.1船測水深數據Fig.4 S&S V19.1 shipborne bathymetry data

對比S&S V19.1模型中616個船測點與筆者收集到的1027個船測點空間分布如圖5(a)所示,圖中“·”號表示筆者收集的船測點位置,“?”號代表S&S V19.1中船測點位置。圖5(a)船測點空間分布比較就總體而言,筆者收集的船測數據與S&S V19.1在大部分區域重合;部分海區S&S V19.1相比于本文增加了多條測線的數據,如研究海區北部,S&S V19.1中包含了筆者不掌握的船測水深數據。需要說明的是,因S&S V19.1空間網格大小為1′,從而基于S&S V19.1提取的船測點間距也是1′;而筆者從NGDC下載的聲學水深測量結果為不規則離散分布,故認為S&S V19.1使用了本文以外的其他水深測量數據,數據量更加豐富。從而為盡可能充分利用可用船測數據,將S&S V19.1中包含的船測點(616個)和本文收集的船測控制點(847個)合并以補償非反演頻段海底地形信息(共1463個)。最終,試驗海區控制點和檢核點空間分布如圖5(b),其中三角形表示檢核點位置,白色六邊形表示控制點的空間位置。

圖5 船測水深點分布Fig.5 Distribution of shipborne bathymetry points

2.2 海底地形模型構建

海底地形主要影響海面重力中長波段信息,因而利用衛星測高重力數據也主要恢復海底地形中長波段信息。依照文獻[31]和筆者前期研究相關結論,本文選擇15 km和160 km作為截止波段,即依托衛星測高重力異常數據主要恢復15～160 km波段海底地形信息,15 km以內和160 km以外的海底地形信息依托船測海深進行補償。首先,利用GMT(generic mapping tools)中格網化操作模塊(surface)將試驗海區離散船測控制點格網化為網格大小為1′的海深格網[32];然后,利用文獻[21]基于撓曲補償和信噪比設計的高通、低通和帶通濾波器對格網海深和S&S V32.1重力異常實施濾波操作,分別獲得長波(波長大于160 km)和短波(波長小于15 km)海底地形,反演波段(15～160 km)重力異常,相應的數據空間分布和統計結果如圖6和表3所示。需要說明的是,依據海底地形解析反演算法[33]可知,重力異常經帶通濾波并向下延拓后與海底地形表現出較強的相關關系,因而為充分利用重力場元與海底地形間相關關系進行統計分析,圖6(c)表示的反演波段重力異常為經帶通濾波并延拓到平均海深面的重力異常結果。

表3 數據統計結果

圖6 不同波段數據Fig.6 Different band data

對比圖6(a)和圖6(b)可發現,長波海底地形表現一種大尺度的低頻特征,而短波海底地形數據則呈現出高頻變化特征,從而可利用重力梯度和重力異常對不同頻段地形敏感程度進行聯合反演[34]。依據式(2)恢復海底地形信息需解決的關鍵技術問題是地形-重力協方差矩陣和重力協方差矩陣的獲取,其中重力異常協方差矩陣為對稱矩陣,地形-重力協方差矩陣是非對稱矩陣。重力異常協方差矩陣求解技術目前已十分成熟[35],而地形-重力協方差矩陣研究成果相對較少。分析地形-重力協方差矩陣元素結構,結合互協方差性質易知,地形-重力協方差矩陣中上三角矩陣元素值可將對應地形和重力位置關系施以地形-重力協方差函數(如本文選擇的似高斯協方差函數)得到;地形-重力協方差矩陣下三角矩陣元素值則是對應重力和地形位置關系施以重力-地形協方差函數獲取。

基于以上分析,首先以反演波段重力異常格網數據為輸入,統計分析試驗海區重力異常分布規律,為使統計結果充分反映統計數據間相互關系,參與統計樣本應盡可能豐富,從而本文實施數據統計分析時,數據間橫向距離和縱向距離最大設置為40′?？紤]到協方差函數非負性質,去掉協方差統計為負值結果,最終得到的重力異常協方差統計結果如圖7(a)所示,距離為零時,協方差統計結果為2 721.248 2 mGal2。依據重力異常協方差統計結果,根據式(11),采用最小二乘方法擬合參數結果如圖7(b)所示,圖中方點為協方差待擬合數值分布,直線為擬合結果,相應的參數計算結果為0.207 8。采用相同的統計策略,得到的重力-地形互協方差和地形-重力互協方差結果如圖8和圖9所示。其中距離為零時,重力-地形互協方差和地形-重力互協方差統計值大小均為33.776 7 mGal·km,其中重力-地形互協方差似高斯函數參數擬合值為0.210 8,地形-重力互協方差似高斯函數參數擬合值為0.192 1。最終,基于似高斯協方差函數的(互)協方差函數表達式為

(12)

圖7 重力異常協方差結果Fig.7 Gravity anomaly covariance results

圖8 重力-地形協方差結果Fig.8 Cross covariance results between gravity and terrain

圖9 地形-重力協方差結果Fig.9 Cross covariance results between terrain and gravity

研究海區格網點間歐氏距離(距離單位為′)解算結果如圖10所示,對圖10距離矩陣施以式(12)重力異常似高斯協方差函數解算,可得相應的重力異常協方差矩陣(如圖11(a))。依前面分析知,地形-重力協方差矩陣為非對稱矩陣,從而將距離矩陣所有元素直接施以地形-重力協方差函數計算結果并非地形-重力協方差矩陣。對于地形-重力協方差矩陣求解,分為兩個步驟:①對圖10距離矩陣的上三角所有元素施以地形-重力似高斯協方差函數,距離矩陣的下三角矩陣所有元素施以重力-地形似高斯協方差函數;②將各自計算結果合并可得地形-重力協方差矩陣,結果如圖11(b)所示?；讷@取的試驗海區地形-重力協方差矩陣和重力異常協方差矩陣,以反演波段重力異常為輸入,采用最小二乘配置方法可恢復相應波段的海底地形,然后疊加依據稀疏船測數據獲取的非反演波段海底地形可獲得試驗區域海深模型。試驗海區海深模型結果如圖12所示,為描述方便,稱為BAT_LSC_1模型。

圖10 距離矩陣Fig.10 Distance matrix

圖12 BAT_LSC_1海深模型Fig.12 BAT_LSC_1 bathymetry model

2.3 模型精度評價分析

引入ETOPO1和DTU18海深模型評估BAT_LSC_1海深模型效能;同時采用未考慮地殼均衡補償的重力異常導納函數理論[36]恢復目標海區海深模型,建立的海深模型稱為BAT_AF_1模型,研究海區各海深模型數值統計結果見表4。

表4 海深模型統計結果

表4中各海深模型統計結果顯示,研究海區海深最深約為5500 m,其中BAT_LSC_1最大海深處海深達到了6000 m;各模型海深平均值為4700 m左右;海深標準差均處于750 m左右,說明本文利用重力異常數據反演海深在試驗海區與目前國際發布的海深模型表現相似。進一步評價本文海深模型質量,以事先準備的外部實測海深數據作為外部檢核參考,采用雙線性插值方法將海深模型值內插到檢核點處,然后與檢核點處海深作差并統計差值結果,各海深模型檢核精度統計結果見表5,表中相對精度為檢核均方差與檢核點平均海深絕對值之比。

表5 海深模型檢核差值統計結果

表5檢核差值統計結果顯示,ETOPO1海深模型檢核差值最大值達到了900 m,其余海深模型檢核差值最值為500 m。就差值標準差和均方差統計指標而言,DTU18在研究海區檢核精度最高,約為百米;BAT_LSC_1和BAT_AF_1檢核精度優于ETOPO1模型,相較于ETOPO1,BAT_LSC_1檢核精度提高了38%左右;BAT_AF_1檢核精度提升約37%;BAT_LSC_1檢核精度與BAT_AF_1相當,檢核精度略高于BAT_AF_1。4種海深模型與檢核點相關系數均達到了0.94以上,其中BAT_LSC_1和BAT_AF_1與檢核點相關系數均高于0.97。ETOPO1檢核相對精度為5.73%,BAT_LSC_1和BAT_AF_1相對精度均優于4%。綜合海深模型檢核差值統計分析結果,驗證了本文方法構建海深模型技術方案可行,即利用統計方法構造地形-重力互協方差函數方法有效。進一步比對各海深模型檢核效果,統計不同檢核差值(差值取絕對值)范圍內檢核點所占比例,結果如圖13所示。

圖13 不同檢核差值比例Fig.13 The ratio of checkpoints in different ranges of difference

對比不同海深模型不同檢核差值范圍內檢核點比例可以看出,DTU18海深模型明顯優于其他3種海深模型。差值小于120 m左右范圍時,ETOPO1模型檢核點數量多于BAT_LSC_1和BAT_AF_1;檢核差值大于120 m左右時,相同差值情況下ETOPO1模型檢核點數量明顯減少,檢核點數量明顯小于BAT_LSC_1和BAT_AF_1;檢核差值范圍大概120～280 m時,BAT_AF_1檢核點數量多于BAT_LSC_1,之后BAT_LSC_1檢核點數量略優于BAT_AF_1。經以上分析可知,試驗海區利用最小二乘配置方法與利用導納函數(未考慮地殼均衡補償)方法構建的海深模型檢核結果大體相當,說明了本文方法的實際可操作性。

2.4 陌生海區海底地形模型構建

據以上研討分析結果可知,采用本文方法可有效構建目標海區海底地形數值模型,根據所建立的海底地形模型可較為清晰地反映目標海區地形地貌特征。最小二乘配置技術恢復海底地形模型關鍵核心是構造地形-重力協方差矩陣和重力異常自協方差矩陣。依統計分析理論,協方差矩陣反映了目標海區相關數據信息分布的內在規律,因此若某一陌生海區數據分布與已知海區數據分布具有相似特性,則可將反映已知海區數據分布內在特性的協方差函數遷移、推廣、應用于陌生海區,進而實現以重力異常為輸入的無船測數據支撐環境下,運用最小二乘配置方法反演海底地形目的。

選擇西太平洋某5°×5°(157°E—162°E,19°N—24°N)海域作為陌生備選海區,以上述1°×1°(155.2°E—156.2°E,22.5°E—23.5°E)試驗海區為已知海區(以下稱為區域1)。將陌生備選海區S&S V32.1重力異常作為已知數據(圖14),相應的數據統計結果見表6,通過比較、評判已知海區和陌生海區重力異常數據相似性程度,搜尋出備選海區中與已知海區數據特征最為相近區域。本文借鑒影像紋理特征匹配相關知識,首先將已知海區和備選海區重力異常作灰度化處理,得到像素值范圍在[0,255]的標準灰度影像,影像結果如圖15所示。

表6 備選海區重力異常

圖14 重力異常Fig.14 Gravity anomaly

圖15 標準灰度影像Fig.15 Standard grayscale image

為確定陌生備選海區中與已知海區(區域1)相似區域準確位置,同時兼顧影像匹配效率,本文分粗匹配和精匹配兩個步驟環節在備選海區完成相似性區域搜索。首先,以已知海區重力異?；叶扔跋窈湍吧鷤溥x海區重力異?；叶扔跋駷檩斎?根據SAD粗匹配算法模型,在備選海區中搜尋與已知海區最為相似的5幅子圖,尋出的子圖位置如圖16中小圖的方框區域。然后,按照灰度共生矩陣(GLCM)和灰度-梯度共生矩陣(GGCM)生成原則,構造已知影像和5幅子圖的GLCM和GGCM,考慮到計算消耗,將灰度值壓縮為16個灰度級(具體構造過程見前文描述),分別選擇GLCM特征值中和平均、和方差、和熵和反差,GGCM特征值中灰度熵、灰度平均和灰度均方差,以及Hu矩(7個不變矩)一共14個特征值作為影像特征向量構成元素;以設計的特征向量為參考,分別計算5幅子圖與已知影像間的歐式距離,與已知影像距離最小的子圖即為最相似海區,精匹配結果如圖16(a)中藍色方框所示區域(以下稱為區域2)。另外,采用以上粗匹配和精匹配方法獲得與已知海域相似度最低(與已知海域歐式距離最遠)區域如圖16(b)中藍色方框所示(稱為區域3)。

圖16 影像匹配結果Fig.16 Image matching results

從S&S V32.1重力異常模型中提取出區域2和區域3衛星測高重力異常數據,相應的統計結果見表7,區域2和區域3重力異常平均值分別為4.13 mGal和25.36 mGal,標準差f分別為43.20 mGal和58.96 mGal。筆者收集到區域2實測水深數據分布如圖17中紅色圓點所示,共計265個;S&S V19.1海深模型中通過實測水深數據約束生成的模型點共309個(認為S&S V19.1中實測水深值約束生成的模型點為實際水深測量值),分布情況如圖17所示,相應的統計結果見表7?？梢钥闯?S&S V19.1中實測水深數據相較于筆者收集的船測數據更為豐富,因此將S&S V19.1實測水深數據作為評估海深模型精度的外部檢核參考。區域3 S&S V19.1實測水深數據統計結果見表7。

表7 區域2、3數據

圖17 S&S V19.1實際船測水深數據和NGDC船測水深數據分布比較Fig.17 Comparison of shipborne bathymetry data between S&S V19.1 and NGDC

據以上研究分析,因區域1和區域2海底地形地貌相似程度高,以區域2反演波段(15～160 km) S&S V32.1重力異常為輸入,將區域1事先解算的地形-重力協方差函數、重力-地形協方差函數和重力異常協方差函數遷移于區域2,應用最小二乘配置技術可獲得區域2反演波段海底地形;而后以S&S V19.1作為區域2先驗海深模型,并施以濾波操作,從而可補充非反演波段(波長小于15 km和波長大于160 km)海底地形信息;疊加區域2所有波段地形信息,進而可實現陌生海域(區域2)海底地形模型構建。采用以上海底地形構建策略,建立的區域2海深模型(稱為BAT_LSC_2)效果如圖18所示。

圖18 BAT_LSC_2模型Fig.18 BAT_LSC_2 model

將區域1相關協方差函數應用于區域2,采用最小二乘配置方法恢復的區域2海底地形模型(BAT_LSC_2模型)可明顯反映區域2的海底地形地貌特征, 區域2海區海底分布著5座大型海底高山。區域2海底地形的呈現狀態定性表明:以海洋表面重力異常為相似性參考,將已知海區解算的協方差函數遷移、推廣并應用于相似未知海區,采用最小二乘配置方法恢復得到的海底地形可較清晰地表達未知海區海底地形特征。為進一步評估該方法的適用性,依據重力異常導納函數反演理論,建立區域2海底地形模型(下文簡稱BAT_AF_2模型)。以事先設置的外部船測檢核點為參考,各海深模型的檢核統計結果見表8。

表8 海深模型檢核結果

表8統計結果表明,BAT_LSC_2檢核統計結果指標優于其余模型,如以檢核均方差為例,BAT_LSC_2檢核均方差為254.07 m,而BAT_AF_2、ETOPO1和DTU18分別為386.62、454.87和308.66 m。BAT_LSC_2與檢核點相關系數最高,達到了0.970 6。4種海深模型中,ETOPO1模型檢核相對精度超過10%,其余模型相對精度均在10%以內,其中BAT_LSC_2模型相對精度為6.16%?；谝陨虾喴治?進一步定量闡明了本文將已知海區相關數據統計計算結果遷移、推廣并應用于相似陌生海區方法可行,構建的海底地形模型具有較強表現能力。比較各海深模型不同檢核差值范圍內檢核點數量占檢核點總數情況,結果如圖19所示。

圖19 不同檢核差值范圍比例Fig.19 The ratio of checkpoints in different ranges of difference

圖19各海深模型不同檢核差值范圍檢核點數量占比曲線比對可以發現,DTU18模型在差值小于400 m左右范圍表現優于BAT_LSC_2,而BAT_LSC_2表現整體優于BAT_AF_2和ETOPO1;BAT_AF_2和ETOPO1模型表現相當。經以上分析可知,就模型檢核統計各指標和不同檢核差值范圍檢核點比例而言,利用最小二乘配置方法構建的區域2海底地形模型明顯好于基于重力異常導納函數構建的海底地形模型,進一步說明了本文解決無船測水深數據支撐環境下海底地形模型構建具有可行性。

類似地,將區域1解算的地形-重力協方差函數、重力-地形協方差函數和重力異常協方差函數遷移于區域3,應用最小二乘配置技術亦可反演獲得陌生區域3海底地形信息,稱為BAT_LSC_3模型。同理,采用導納函數方法獲得的區域3海底地形模型命名為BAT_AF_3模型。BAT_LSC_3模型和BAT_AF_3模型統計結果及其外部水深檢核統計結果分別見表9和表10。

表9 海深模型

表10 海深模型檢核結果

表9顯示BAT_LSC_3模型和BAT_AF_3模型整體表現基本一致。由表10可知,BAT_LSC_3模型和BAT_AF_3模型檢核均方差分別為347.03 m和359.56 m,BAT_LSC_3模型精度略優于BAT_AF_3模型。相較于BAT_LSC_2,BAT_LSC_3質量明顯下降:①BAT_LSC_3系統誤差明顯增大,BAT_LSC_2檢核差值平均值為46.08 m,BAT_LSC_3檢核差值平均值達到了182.73 m;②BAT_LSC_2檢核精度明顯優于導納函數結果,而BAT_LSC_3與導納函數海底地形結果比較接近。不同檢核差值范圍內BAT_LSC_3模型和BAT_AF_3模型檢核點數量占檢核點總數情況如圖20所示。

圖20 不同檢核差值范圍比例Fig.20 The ratio of checkpoints in different ranges of difference

對比圖20和圖19可知,BAT_LSC_3模型與導納函數反演結果較為接近,檢核差值小于250 m左右,BAT_AF_3模型表現明顯優于BAT_LSC_3;而BAT_AF_2模型則在整個檢核差值范圍表現均好于導納函數結果。筆者認為原因是:使用最小二乘配置方法構建區域2和區域3海底地形模型時,均使用了區域1地形-重力協方差函數、重力-地形協方差函數和重力異常協方差函數進行解算,而區域2和區域3分別為與區域1相似度最高和最低區域,從而區域2協方差函數與區域1更為接近,進而區域2解算結果質量更高。另外值得說明的是,本文利用重力異?；謴秃５椎匦蔚淖钚〗刂共ǘ螢?5 km,即奈奎斯特(Nyquist)波長為15 km;則依據頻率域采樣定理可知,恢復的海底地形分辨率最高為30 km,而本文構建的海底地形模型在波長小于15 km部分利用了S&S V32.1模型補充,因此構建的海底地形模型最終分辨率依賴S&S V32.1模型,海底地形網格大小為1.8 km左右。

3 結 論

本文嘗試探索基于已知海區數據統計信息恢復未知陌生海區海底地形的手段途徑。首先,利用已知海區可用的稀疏船測水深數據和衛星測高重力異常數據,統計建立了數據已知海區(區域1)地形-重力互協方差函數和重力自協方差函數,依據最小二乘配置方法恢復了區域1的海底地形;然后,以衛星測高海面重力異常作為相似性參考,采用影像粗匹配和精匹配算法模型,尋出了與已知海區重力異常最為相似的無船測水深數據海區(區域2);最后,將區域1相關數據的協方差函數遷移、推廣、應用于區域2,采用最小二乘配置方法構建了相應的海底地形模型,并評估了所構建的海底地形模型效能。經過定性和定量兩個層面比對分析試驗結果可以發現:本文方法反演的海底地形模型可有效描述研究海區海底地形地貌特點,反演的研究海區海深模型整體呈現效果好于根據重力異常導納函數構建的海深模型,同時本文方法可巧妙回避導納函數涉及的地球物理參數獲取難題。試驗區域1中,BAT_LSC_1模型檢核精度相較于ETOPO1提升了約38%;試驗區域2中,無論檢核點差值統計結果,還是不同檢核差值范圍檢核點比例情況,BAT_LSC_2表現均明顯優于BAT_AF_2和ETOPO1。